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Capacity of the Gaussian Channel With 
Memory: The Multivariate Case 


By L. H. BRANDENBURG and A. D. WYNER 
(Manuscript received November 13, 1973) 


A formula ts derwed for the capacity of a multi-input, multi-output 
linear channel with memory, and with additive Gaussian noise. The 
formula 1s justified by a coding theorem and converse. The channel 
model under consideration can represent multipair telephone cable in- 
cluding the effect of far-end crosstalk. For such cable under large signal-to- 
noise conditions, we show that channel capacity and cable length are 
linearly related; for small signal-to-noise ratio, capacity and length are 
logarithmically related. Crosstalk tends to reduce the dependence of 
capacity on cable length. Moreover, for any channel to which our capacity 
formula applies, and for large signal-to-noise ratio, there 1s an asymp- 
totic linear relation between capacity and signal-to-noise ratio with 
slope independent of the channel transfer function. For small signal-to- 
noise ratio, capacity and signal-to-noise ratio are logarithmically related. 
Also provided 1s a numerical evaluation of the channel capacity formula, 
using measured parameters obtained from an experimental cable. 


l. INTRODUCTION AND STATEMENT OF RESULTS 


Our problem is to calculate the capacity of a multi-input, multi- 
output linear channel with additive Gaussian noise, and to justify the 
formula by a coding theorem and converse. Specifically, we consider 
the following channel. The channel input and output are sequences 
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{fa(n)} 22, {y(n)}2. of real s-vectors* related by 
y(n) = h(n — be(k) + a(n), (1) 


where {h(m)}n=—0 is a fixed sequence of real s X s matrices (the 
indicated operations being ordinary matrix arithmetic), and {z(n)}2. 
is a sequence of Gaussian random s-vectors for which Ez(n) = 0, and 


Ke(n)Le(n as m) }* =n r(m), —-eo<cn m< m, (2) 


where 7(m) is an s X s matrix. The motivation for this problem is that 
it is a model for a multipair telephone cable. 

The first sections of this paper are highly theoretical; the formula 
for channel capacity is carefully and precisely established by means of 
several rather technical theorems. In the final section, Section IV, we 
discuss some engineering implications of our formula in terms of its 
asymptotic behavior, and evaluate the capacity numerically with 
measured parameters obtained from an experimental multipair 
telephone cable. 

A code for this channel with parameters (M, N, S, d) is a set of 
pairs { (Xi, B,) han, where xX; = ( es) x(—2), x= 1), x;(0), a;(1), _ aS 
is a sequence of s-vectors that satisfy the following: 


zi(n) = 0, for n<0, and n2QN, (3a) 
| N=1 
+ lledn)ll? < 8, (3b) 
n=0 
(where ||-|| denotes Euclidean norm) and the B; are (measurable) 


subsets of ®*% with the following property. Let y; = (---, yj(—1), 
yi(0), yi(1), ---)* be the channel output vector which results when the 
channel input is X;, 1.e., 


yi(n) = > h(n — k)ae(k) + 2(n) 


N-1 
Pa h(n — k)ax:(k) + 2(n). 
Let yi” = [yi(0), ---, y(N — 1)]' © a. Then Bi(1 S$ 7 S M) must 
satisfy 
Pe & Pry © Bi} Sd. (4) 


*Vectors will be taken to be column matrices unless otherwise indicated. 
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Thus the x; are code words and B;is the set of output N-vectors 
which are decoded at the channel output as x;. Inequality (4) expresses 
the requirement that the error probability, given that x; is transmitted, 
does not exceed A. 

A number p 2 O is said to be ‘‘S-admissible’ (S 2 0) if for all 
\ > 0 there is an N such that there exists a code with parameters 
([2%-], N, S, A). The channel capacity C's is defined as the supremum 
of S-admissible rates. Our problem is the calculation of Cs, which we 
shall solve provided the channel satisfies the following conditions: 


(2) The filter {h(m)}: We assume that the filter is causal, i.e., 
h(m) = 0, m < 0. We also assume that 


2X [html < ©, (5a) 
and that there exists a B > 0 such that for m > 0, 
[h(m)|| S$ Bm-, (5b) 


where the Euclidean norm ‘“‘||-||’’ of a matrix is the square root 
of the sum of the squares of its entries. From (5), the (discrete) 
transfer function, 


H(@)= > bine, —2 S67, (6) 


exists and is continuous. H(@) is an s X s matrix. We assume 
that for —7 S 0 S a, det H(0) ¥ 0. 

(11) The notse covariance: We assume that the covariance sequence 
r(-) satisfies 


X Ir@ll < @, (7) 


so that the (discrete) power spectral density 


ITA 


R(0) = >> r(nye®, —2 < 6 


neem 


T (8) 


exists and is continuous. R(@) is an s X s matrix. We also 
assume that for —7r S @ S zw, det R(6) € O. 


We can now give the capacity formula. Let the s X s matrix! 
T(6) = H(@)-1R(6)H(@)-*, and let \1(@), \2(@), «++, Ae(6) be the eigen- 


t For any nonsingular complex matrix A, A~* is the transpose conjugate of A7}. 
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values of (0), -7 S 6S mw. Then),;(¢) >0,1 Sj Ss,-7rSO08 7. 


Let S = 0 be given, and let Kg be the (unique) positive number such 
that 


= [* domax[0, Ks — (6)] = 8. (9a) 
a ee ee 
Then 
Cs = Eas , d6 max (0 logs 5). (9b) 
Agr j=l —7 ; d; (8) 


Our main result is the following. Consider the channel defined by 
(1) with H = H,(@) and R&R = F,(@). Then Cg is calculated from (9) 
with I'(@) = H,(0)—R1(6)H1(6)-*. 


Theorem 1a (Converse): Let p= O be an S-admissible rate for this 
channel. Then p S Cs. 


Theorem 1b (Direct-Half): Let S20, «€> 0 and p(O Sp < Cs) be 
arbitrary. Then for N sufficiently large, there exists a code with parame- 
ters (M, N, S, d) where 


Mz2=ee% and \S.«. 


Sections II and III of this paper are concerned with the proof of 
Theorem 1. Section IV is concerned with the asymptotic behavior and 
numerical evaluation of the channel capacity formula (9), with specific 
attention to multipair telephone cable. 

Theorem 1 is very similar to the results on continuous-time Gaussian 
channels due to Holsinger and Gallager.! In fact, for the special case 
in which H(6) is the s X s identity, the theorem follows immediately 
from the analysis in Ref. 1. We suspect that it might be possible to 
obtain all of Theorem 1 by paralleling Gallager’s techniques for this 
discrete case, although such an approach is somewhat more cumbersome 
than the approach followed here. Furthermore, the present approach 
lends itself immediately to broadening the model to consider the effects 
of intersymbol interference from previous channel uses, as Gallager’s 
approach does not.” In fact, to establish Theorem 1b for the inter- 
symbol interference channel we require only to add in one of our 
lemmas a term ‘yy; (representing the effect of previous channel 
uses), and to show that its norm |||¥ys3|/| = o(2V?). 

Careful analysis of the proof of Theorem 1 will indicate that the 
conditions on the filter {h(m)} given in Section I can be replaced by 
simply requiring that the filter have a causal inverse {g(m)}, such that 
g(m) = 0, m = mp (i.e., finite memory). Thus, our results contain a 
generalization of those given by Toms and Berger.’ 
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Il. NOTATION AND MATHEMATICAL PRELIMINARIES 


Let I§(a,b), s=1,2,---, —- 7 Sa<bS ~, where a, b are 
integers, be the set of sequences {x(n)}4.4 where x(n) is a real s-vector. 
Such sequences will often be written as column matrices x = [2*(a), 
x(a +1), ---, 2*(b) ]*. Let ||-|| denote ordinary Euclidean norm in 
s-space, and for s X s matrices A, let ||A|| be the Euclidean norm (i.e., 
the square root of the sum of the squares of the components of A). 
For sequences x € /$(a, b), the (Euclidean) norm is 


lll = [2 eon | (10) 


The space /§(a, b) is a Hilbert space with the obvious inner product, 
written (x, y) = Din=a2'(n)y(n), xy € 1"(a, 6). For x € lp (— ~, o), 
denote by x“) € 1§°(0, N — 1) the column matrix 


x1) = [x‘(0), x'(1), ng a'(N = eae (11) 


We denote operators on I$? (a, b) by script letters, e.g., 5. We define 
the norm of § by 
5) = sup Sz 


12 
reine) [IFT — 





If |F| < ©, we say that F is bounded. An operator § is said to be of 
a convolution type if, for x € l§(a, b), 


b 
(Fx) (n) ane py f(n a k)x(k), maar"; b, (13) 
where { f(n)}2=3 is a fixed sequence of s X s matrices and the indicated 


operations are ordinary matrix arithmetic. Let L be the set of convolu- 
tion-type operators on 1§?(— «©, «) for which 


y Isl < ©, (14) 
For operators § in L the transfer matrix 
F(o=) = > fine, -2 S057 (15) 


is well defined. F(6) is ans X s matrix and is continuous (in Euclidean 
norm) for —r S@X 7. Concatenation of operators 51, S2€ L, 
defined by sequences { f:(n)} and { f2(n)}, results in a convolution-type 
operator 53 = S2:5, in L defined by the sequence { fs(n)} where 
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fa(n) = Diva—w fo(k) fi(n — k). Further, the corresponding transfer 
functions satisfy /3(0) = F.(@)F1(@). Other relevant properties of the 
class L are given in the theorem below, the assertions of which are 
generalizations of well-known scalar results. The proof is given in the 
appendix (Section A.1). 


Theorem 2: Given § € L, defined by { f(n)} or F(@), then 
a. || 
b. || 


IA 


Dn=-« ||f(n)|| < ©, so that F is bounded on If) (— 0, &), 


IIA 


_max |[F'(9)I 
c. If & ts self-adjoint, then for all x € IS (— ~, ~), 


E at(n) s(n — m)(m)| = (_max_|PCII)Ilxil 


d. § has a bounded inverse denoted 5-1 if and only rf det F(0) ¥ 0, 
—a S087, in which case the transfer matrix corresponding to 
5-1 7s [F(6) |-1, and $ € L. 


Let z = [---, 24(—1), 24(0), 2(1), ---, ]* be a sequence of random 
s-vectors with covariance 


Ez(n)z(n — m)* = r(m), 


where r(m) is an s X s matrix. Under the assumption that 


x lir(m)|l < @, (16) 
the power spectrum 


R(@) = dX r(mem, —2SO0K7, (17) 
is well defined. Let 5 be an operator in L corresponding to the transfer 
matrix F(6). ThenzZ = $zis asequence of random s-vectors with covari- 
ance {7(m)} and corresponding power spectrum R(@) = F(6) R(6)F(6)*. 

Let z be a sequence of zero mean Gaussian n-vectors with covariance 
r(m) satisfying (16) and power spectrum R(@). Let Amin(@) be the mini- 
mum eigenvalue of the matrix R(6). Let z“ = [2#(0), 21), °°, 
z' (N — 1) ]'be a segment of z of duration N. Then there exists an 
(N-s) X (N-s) matrix Ty such that 


w= Tyz™), (18) 


is “white,” i., Hww'! = Iy.,. The indicated operation in (18) is 
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matrix multiplication. The only property of Ty which we need here 
is that for any NV -s-vector u 


Poulll s | ey | lull (19) 


—7rSO0S9r 


Finally, let 3¢ be the convolution-type operator on 1§)(— «0, «) 
defined by {h(m)} (Section I). Note that by (5) 3 € L, so that by the 
assumption following (6) and by Theorem 2d it has an inverse, say 
G = 3c~, of the convolution type in L. Let {g(n)}". be the sequence 
which defines G. Let the operator Gy(N = 0, 1, 2, ---) be defined by 


(Gex)(n)= F g(n— Balk), —e<n<e. (20) 


We conclude this section by stating as lemmas two known results. 
We explain in the appendix (Section A.2) how to obtain these results 
from published material. 


Lemma 8: For the special case when H(@) = I, (the s X s identity) and 
R(@) = R2(6) [ so that T = H-RH-* = R2(6) ], say that x ts a random 
channel input sequence for which x(n) = 0,n <0, 2 N, and H]\|\x{||/? 
< NS(S>0, N = 0, 1, 2, ---). Let y be the corresponding channel 
output sequence. Then, the mutual information 


I{x,y} SNCs 
(where C's ts calculated with T(@) = R2(@)). 


Lemma 4: For the special case where H(0) = I, and R(@) = R2(6), then 
we have a stronger version of Theorem 1b: Let S = Oand p (0 S p < Cs) 
be arbitrary, where C'g is calculated with T(@) = R2(0). Then, for N=1, 
2, ---, there exists a code with parameters (7, N, S, X) where 


M2=2e" and XS Ae®*%, A, B>O. 


Ill. PROOF OF THEOREM 1 
3.1 Converse 

Let {(x;, B;)}7’ be a code with parameters (7, N, S, A) for the 
channel of (1) with H = A,(@) and Rk = R,(6). Let x be the random 
sequence which results when x = x; with probability 1/M (1 S27 MM). 
Let y = 3x +z be the corresponding output sequence, and y‘¥? 
= (y(0)', ---,y(N — 1)*)*. The theorem will follow in the standard way 
from the Fano inequality (see, for example, Ref. 1) if we can show that 


I{x, y} < NCs. (21) 


GAUSSIAN CHANNEL WITH MEMORY 751 


But 
E{x,y} s I{x,y} = I{x, Hy} = If{x, x + 2}, (22) 


where Zz = iz is a stationary Gaussian random process with power 
spectrum I'(@) = H;'(¢)R.i(6)H,~*(6). Thus, we can apply Lemma 3 
(since x satisfies the required hypotheses) with R2(@) = I'(@) to ob- 
tain (21). Hence, the theorem follows. 


3.2 Direct-half 


Consider again the special case of our channel where H(@) = J, and 
R(6) = R.(@). The idea behind the proof is to construct codes for the 
general case (H, FR arbitrary) by modifying codes (whose existence is 
guaranteed by Lemma 4) which are known to be good for this special 
case. We proceed as follows. Let {(x;, B;)}7* be a code for this channel 
with parameters N = N,; and S = S;. Then, for 1 S 7S M, we have 


yi = x 4 70), 


where the superscript operation is defined by (11). Let Tn be the 
whitening filter (discussed in Section II) for which Tyz‘*) = w and 
Eww't = I[y,;. Letting v; = Tyy{™ and u; = Tyx{”, we have 


Vv; = Uu; + W. (23) 


Let us assume that the {B;} correspond to the minimum distance 
decoder, ie., y” € B; if ||lv— ual|| < ||| v — u,||| for all 7 4 7, 
where v = TJ'yy‘?. Then 


Pec = Priyl? & Bs} = Pr Y lllve — wall [lve — wll 


Pr UY Aillwill > [Ilw — Gas — wdHID 


Ee YU {(w, u; — uy) 2 9||[u; — uill|?}. (24) 


Thus, in particular, for all 7 + 2, 
Poi 2 Pr{(w,u; — us) 2 3||lu; — uyl||?} 
= 8.(3||luy — uill[) = B-L3[||Tw(a — xi)|[[], (25) 


where 


®.(u) = Ao? era 


V2Qa Ju 


is the complementary error function. 
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Let H,(6) and R,(@) be arbitrary, and suppose that we are given a 
code {(x;, B;)}} with parameters (M, N, S, \1) for use on the special 
channel with H = I, and R(é) = R. = I(6) = Ay'RiH7*. Assume 
that the 8; corresponds to the minimum distance decoder so that P.; 
is given by (24). We now construct a new code {(x;, Bi)}*4, with 
parameters N = Ne = (1+ 6)Ni and S = S2 = a’S,/(1 + 8) for 
use on the general channel with H,(6), R1(@) arbitrary. We set 


az(n), OSnS N, —1, 


0, otherwise, 8) 


vi(n) = 
where a > 1 and 6 > O are arbitrary. Note that we have allowed a 
guard-band or dead-space or width 6N, following the channel input 
signal. The decoding sets B; (1 S 7 S M) are described below. 
The channel output is as in (1) 


y = Hx + Z, 


where x is the channel input, z is the noise, and JC is the operator corre- 
sponding to H,;(6). Let Gy, be the operator defined in (20), and let 


y= Gn = Gn,dtx + Gn, 
=x+Z+%& + &, 


where & = Gy,ix — x, Z = Gz, and & = Gy,z — Z. Let us note that 
y is calculable from y‘%») = (y*(O), ---, y(N2 — 1))*. Further, the 
noise Z has power spectrum I'(@) = Hy71(6)Ri(0)H; (0). (In fact, if 
—, = & = 0, the channel would be equivalent to the special case, and 
the direct-half of the coding theorem would follow from Lemma 4. 
Although this is not the case, of course, we will show that &, and £2 are 
sufficiently small so that Lemma 4 can be applied anyway.) The 
decoding sets B; are defined by: y4” € Brif FY) CB, 1 Sis M. 

Letting yi = x; +Z2+8+%& (1 SiS M), and letting Ty, be as 
above, we define 


vi S Ty yi = Ty, xt + Ty) + Ty Be 4 Tye"? 
au; + w+ yi + ¥2, (27) 


where u, and w are exactly as in (28) and y; = Ty,&{""(¢ = 1, 2). The 
decoder for the derived code is the minimum distance decoder for the 
v*’s. Now, following the same steps as in (24), we have 


Pa = Pry & Bi = Pr YU {Illve — anaslll 2 lve — ull 


= Pr {wt va + ya ay — us) 2 5 Illy — wall. (28) 
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Now, according to (9), the channel of Section I [with arbitrary 
H,(@) and (0) |] and the special channel with H(@) = J, and R(@) 
= H,(6)—1R,(6)H,(6)-* = 1'(@) have the same C's. Let e, S > 0 and 
o (0 S p < Cs) given, and let {(x;, B,;)}{ be a code with parameters 
(M, Ni, Si, \1) (as guaranteed by Lemma 4) such that 


M = 2°": and hy Se. 


We will show that with N, sufficiently large, the derived code has 
parameter \ S \i + e. Thus, we will have found a set of codes with 
parameters (7, Ne, Se, \) with 


2 i 
vem a’; M 2 exp] Prat 





(1+ 6)’ 1+6 
and 
Ww Ze, 


Since Cg is continuous in its arguments, and a may be chosen arbitrarily 
close to 1, and 6 arbitrarily close to 0, the direct-half of the coding 
theorem (Theorem 1b) will have been established. 

To show that \ for the derived code S \; + e, we must show that 
for each 2 = 1, 2, ---, 


Pr{yi®® & Bi} S PriyY? EB ty} +e (29) 


Inequality (29) will follow directly from the following lemmas, the 
proofs of which are given at the end of this section. 


Lemma 8: Inequality (29) ts satisfied af 
(a =) 
2 





Pr{{{ly1 + alll S min (liu; — uill[} 21—¢ (30) 


Lemma 6: For the codes {(x;, B,;)}7‘, as N- ~, 


min ||ju; — u,|[/? 2 OC). 
Aj 


Lemma 7: For arbitrary a > 0, 
Pr{\ily: + yall? S aNi} 1, as Nim oO, 


Now, from Lemmas 6 and 7, condition (80) in Lemma 4 will be 
satisfied for N, sufficiently large. This establishes Theorem 1b. 


Proof of Lemma 6: Let 


IIA 


8 = [ills + rll s SS 


By hypothesis, Pr{S} = 1 — .. Since B; and B; correspond to the 


min |||u; — u,l||}- (31) 
VF) 
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minimum distance decoder, we have from (28) 


Pry Oe Bye PS 1) ye EB) Pris 
< Pr Us C) (w a ae Y2, Uj; — U; = . I|{u; = wl + €. (32) 
j rt 


Now if S occurs, 


yi + v2, Us — Ua)| S [Ilva + velll-{[fay — ull] 


— ] 
s (* > )ilius — wile 
Thus, the event in the right member of (32) satisfies 





(w, uy — ws) & Say — alll? — (2+ )iias — wal 


€ {(w, uy — us) 2 3|fuy — ail ||}, 
and (32) becomes 
Priyi®® & BY} S Pr U (dw, uy — ws) & filly — willl) + « 
= Prly € Bi} + ¢ 
where the last equality follows from (24). This is (29) so that we have 
proved Lemma 5. 
Proof of Lemma 6: For the codes { (x;, B,)}{’, and Ni— &, 


P.; Ss es NA, 





8A) fw + rat toy wy — wy) 2 $ lls — wll 
| 


so that from (25), 
&.(3||[u; — uyl|]) S Aew?%. 


Since, as 7 > ©, ®,(7) = e7 (7/2 [+e(D], we have 
[fas — usl[[? 2 8BNiL1 + o (1)], 

which implies Lemma 6. 

Proof of Lemma 7: First note that by (19) 

[Ilya + vrell| 


= [Psa + UL S | ey lige + EFI 


min Am 
an eee | ciiseu + [Nem 
= | min Amin(6) : " 
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o(N}) as Ni > ~, and 
aN,} —0, as Ni— ©. 


I 


V 


Thus, it will suffice to show first, |||E{"”]]| 
second, for arbitrary a > 0, Pr{|{|E&%|||? = 
1. Let x be one of the code vectors in { (x;,B7)}. Then =; = Gy,3Cx—x, 
so that for —~ <n < a, 
ta(n) =X g(n — b) (sex) (h) (33) 


k2=N2 

Now, since x is one of the code vectors {x;}”, x(n) = 0 for 

nELO, Ni — 1]. Also since 3 is causal, we have (3x)(k) = 0, 
k < 0, and 


s(n) = g(m — k)(sex)(R) 
=F gin —b) YE hk— Heals). (8H 


Next, define the sequence t) by 
k= Ne 


alan . . 
Wk) = 1 a hk — g)x(9), 
k< Ne 


Then (34) is 
f(n) =D g(n — ky(h), 
1e., £1 = Su, and 

| Ell] S 1S] fell. (35) 


Now |S| S$ /iw__.q |lg(n)|| < © from Theorem 2d, and 


wlll? = Wweole = 3 |S ace — poco |f 
SCE le - ae\(CE lew), 


< 


where in the final step we have used the following form of the 
Schwarz inequality: if a is a sequence of s X s matrices, and x is 


a sequence of s-vectors, then 


IX a(me(n) | SE llacapll? D le) 
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But since 25? |]z2(7) |]? = |||xl|]? S a2S1Ni, we have 
co «=—sdsNy—-1 
PUT? S atSiW1 2 2 [lace — 9D! (36) 
We will now show that, as VN; > ~, 
co Ny-1 o Ny-1 
2 de WE Pe De a 2, (37) 
k=N2 j=0 k=Nq j=0 


where f(k) = ||h(k)||. Expressions (35), (36), and (37) together 
imply that 
Ex ||? Ss [[fEl[[? = oN), as Nim &, 


which is what we set out to establish. It remains to establish (37). 
Now, from (5), is f(k) < », and f(k) S B/k. Setting 


F(R) = © P), 


we have 
A co N,-1 
Q= 2) Lb KPk-)) 
k=Nq j=0 


k=N2 t=k-N1+1 


p@= [Pe — Ni +1) — FQ) 


N2 N2 
> Fe Sy FO). 


k=No-—Ny+1 1+1 
Now* 
N2 N2 N2 
» Fk) =kF(K)) + De (k-DPKR— 1). 
6Nit+1 5Ny 6Ni+1 
But 


kF(k) i < N.F(6Ni1) = Ne > fk) Ss Gro > kf?(k), 


* We have made use of the formula (summation by parts) 


b 
TM eaua <2 eS ew), 


a—l 


where Au(k) = u(k) — u(k — 1). Here v(k) = F(k) and u(k) = k. 
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758 


and 


2% = DPb-1 s FW, 


5N4 


Thus, 


as (442 +1) Farm s (AE )ay say 0, 

5N1 
since >-¢ f(k) < ». Thus, (37) is established and we have 
finished the first part. 


2. Since for anya > 0, 


< Elle i? 
aN; 


Pr{||[Es" |||? 2 aNa} 
and since Ne = (1 + 6)N,, it suffices to show that 
El\|&3"? |||? = o(N2), as Ne, 


Now & = Gy,z — Z, where Z = 3C—!z. Hence, 


EL Eo(n) &(n) ] = i”) Georg i = 7). 


1j 2N2 


Let 8 denote any fixed s-vector and define a sequence w of s- 
vectors by 


y(n — 1) = 
Then 


gi(n—168, t+<O and +2 Nz. 
0, 0S1< Ne, 
6 BL elm) B(m) B= Lon — Dr — Hvln — 3). 


Since r(-) is a covariance, r(k) = rt(—k) for all k. Application of 
Theorem 2c shows that the double summation above is bounded 
by 

_max ||R(6)|] 2 [v(m — aI). 


Since 


N 3-1 No-~1_ s 
(ee = BOE) = LL ekalnyesn)e, 
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where e, is the s-vector with jth entry 6,; (», 7 = 1, ---, 9s), the 
desired result will follow if we show that 


Ni— 


1 g—1 
a8 %, HL te(n)#4(n)I¢ = 0(1) 
v4 n=0 
for any 8. Thus, it suffices to show that 
]. Nae! 00 ; 
wy, DI — al? = of). 
2 n=0 i=—»o 


But from the definition of y, the last expression can be rewritten 
as 


] Ne-1 00 
wey, &, CV@IP+Io(-B IPD. 
2 n=0 k=n+1 


Since G is in L, Theorem 2a implies wu is in 1 (—«, «). In 
particular, asn— © 


&, Civ@oll?+llv(—B7] = of) 
and the desired conclusion follows immediately. 


IV. ASYMPTOTIC BEHAVIOR AND NUMERICAL EVALUATION OF THE 
CAPACITY FORMULA 

In this final section, we discuss some implications of the channel 
capacity formula given in (9). As in the prior sections, the channel is 
assumed to be a multi-input, multi-output channel with memory, 
and with additive Gaussian noise. But in contrast to the previous case, 
instead of a discrete time channel, we consider an equivalent continuous 
time bandlimited channel. 

Specifically, the channel inputs or code words (in a 7 second block 
coding interval) are vector-valued functions z(-) of dimensions s, 
bandlimited in frequency interval [—W, W ] such that the samples of 
x(-) satisfy 


nN nr (1) 
es = a Sos > 


We also have the average power constraint 
Jo llewllat s 87, 
where ||-|| denotes Euclidean norm. 
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The channel has s inputs and s outputs and has transfer function 
matrix H(f) for fe [—W, W]. The additive noise is also vector- 
valued and has power spectral density matrix R(f). Let {,(f)} 
denote the set of eigenvalues of H-3(f)R(f)[H-"(f) ]*. 

The capacity of this channel is determined as follows. Let S > 0 
be given and let K be the unique positive number satisfying 


ah [- max[0, K — \,(f) Jdf = S. (38a) 
Then 
12 Ww K 
C= 5 a [. max( 0, loge KP) Jas (38b) 


with C in bits per second. 

Formula (38) can be obtained from the analogous formula (9) via 
application of the sampling theorem. The somewhat tedious derivation 
is carried out for the scalar case (s = 1) in the appendix (Section A.3). 

We consider several implications of (88). Specifically, for large 
signal-to-noise ratio, C' is linearly related to signal-to-noise ratio; a 
change in C is proportional to a change in signal-to-noise ratio (in dB). 
Furthermore, the constant of proportionality depends only on the 
product sW and is independent of any other characteristic of the 
channel. For small signal-to-noise ratio, Cis logarithmically related to 
signal-to-noise ratio; a change in logioC' is proportional to a change in 
signal-to-noise ratio (in dB). The constant of proportionality is 0.1 
for any channel. In the case in which the channel represents multi- 
pair telephone cable with small far-end crosstalk, we show that for 
large signal-to-noise ratio, C is linearly related to the length of the 
cable, and for small signal-to-noise ratio, C is logarithmically related 
to length. Furthermore, the effect of the crosstalk is to reduce the 
dependence of C on cable length. Finally, we present a numerical 
evaluation of (38) using realistic parameters obtained from an experi- 
mental cable consisting of two twisted pairs of wire. 


4.1 Dependence of channel capacity on signal-to-noise ratio 


Then N, represents the noise power per hertz, per dimension, and 
sWN, represents the total noise power. We define the following 
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normalized quantities: 
Nf) = 2A,(f)/No 
K* = 2K/N. 
P= S8S/sWN, 
P* = 10 (logio P). 
Now P* is a measure of signal-to-noise ratio in dB. By substituting 


_ the above quantities into (38), and using the fact that each ); is a 
symmetric function, we have 


P=— = * ea: k* — vi(f))df (39a) 


and 


s Ww 

G2 (0, 7 39b 

» 0 ees ee Ni rif) J ( ) 

We will determine the asymptotic behavior of C for both very large 

and very small P. For this purpose we define for every number K* 
sets A;,7 = 1, ---s as 


= {f:\(f) S K*; f = 0}. 
Let 6; be the measure of A; and define 6 = (1/s) >> 6;. In addition, we 
require the definition of two average channel characteristics, \ and 
log X. Let 


>| 
ll 


j ee * 
5, |, MAS, 
and 


a, Le , 
logh= = [log ri Nas. 
SO{=1 A; 


Note that 6, \ and log ) are all functions of P. Let Amin = min{d;(f): 
0s fs W;1 8218 s}. Recall from Section I that \nin > 0. 
Now, from (39), 


and 


= loge K* — log X. 


These equations combine to yield 


C 


WP ee 
= loge (= + 1) + loge \ — log X. (40a) 
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We investigate (40a) for large P. Assume all the \7’s are bounded. 
(Actually, boundedness follows from the hypotheses in Section I.) 
From (39a), K* is an increasing function of P. Let P be sufficiently 
large so that \,(f) S AK* for all f€ [0, W] and 7 = 1, ---s. Then 
5 = W and (40a) yields 

C P | aes 

We logs (5 3 1) + (loge A — log A). (40b) 
Note that for given i, s, W, and P, C is minimized when all the };’s 
are equal and constant. Now, for P> A, we have from (40b), 


CX sW (loge P — log j), 


or 
CH 0.3322 sWP* — sW log X. (41) 


Now (41) represents a line in the C — P* plane with intercept 
—sW log d and slope 0.3322 sW, and the region of validity of (41) is 
P> . Note that the slope is independent of the \7’s; the intercept 
and region of validity are determined by the average channel character- 
istics \ and log \ evaluated over the whole interval [0, W]. 

We now investigate (40a) for small P. Observe that (WP/5\) + 1 
= K*/), and as P approaches zero, both K* and ) approach Amin. 
Hence, WP/6\ — 0, as P > 0. Then (40a) is approximately 


CWP ._ 
5h gy 1082 & + loge d — log a, 


which can be rewritten as 


(loge A — log d) ce 


C 
7 © | lows e + Keay x 


Ww (40c) 


We show in appendix Section A.4 that for any channel characteristic 


NN al 


. logedA — log dX 

] a a IT 0. 42 

pao K*¥— 
Hence, for small P, 

(loge e)sWP 
CR ST, 
Amin 
or in logarithmic terms, 
P* 
logis C 10 + logis (SW loge e) — logio Amin. (43) 
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Now (48) represents a line in the logio C — P* plane with intercept 
logio (SW loge e) — logic Amin, and slope 1/10. However, the region of 
validity of (48) is difficult to specify because the location of this 
region depends not just on Amin but on the shape of the channel charac- 
teristic in a neighborhood of Xin. 


4.2 Dependence of channel capacity on cable length 


Suppose that the channel characteristic is a function of a length 
parameter J. Let /; and J, be two values of J and suppose that P* is 
large and the channel capacity vs P* characteristic is in the linear 
region for J; and J». Then, from (41), 


C(lz) — Chi) & sWLlog Ai) — log A(l2) J, 


or 





ee " Ni(la; f) | 
C(l) — CL) Ef lows (Se Jas (44) 


where in these relations we have explicitly shown the dependence of 
\; on length. If P* is very small so that (43) is valid, then 





Nmin l 
leew) = locus = logo( Amin (U1) ). (45) 
Amin (C2) 


Now consider a multipair cable of length J, with s twisted pairs, 
small far-end crosstalk, and additive white noise. We assume that the 
crosstalk voltage on a single pair due to all disturbers is proportional 
to l? f. Assume also that the attenuation on any pair is proportional to 
lf?. If the crosstalk is very small, then a reasonable form for }; is 

ebilf? 


ML; f) = Tf edf? 


where 6; and c; are constants related to attenuation and crosstalk 
coupling.* Define the averages b and cas b = 3° b,/s and c = >. c/s. 
Now ; can be expressed as 

Mi = exp[b,lf? — In(l + eo f7l)], 
and for small crosstalk, we have c;f?] « 1 for all 7 and all f and / ina 
range of interest. Then we have approximately 

rn; ~ el(bificif?), 

and, from (44), 


AC Sf Dg : 
ar ~ towel) & fo if - of )ds, 
~=1 J/0 
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which can be evaluated as 


AC ‘ cWi 
Al 0.96 bsW (1 r ) (46) 


(for large P). Thus, C and 1 are linearly related. Note that the effect 
of the crosstalk is to reduce AC’/Al. If cW?/2b & 1, then C is effectively 
independent of length. It is theoretically possible to have cW?/2b % 1 
and yet have cWl «<1 as required by our analysis. These relations 
imply that 2W#lb « 1 is a necessary condition that very small cross- 
talk significantly reduce AC/Al. However, we expect that for realistic 
cable parameters, the reduction in AC/Al due to small crosstalk will 
not be significant. 

Now assume that the channel does not pass dc; i.e., the channel 
characteristic is that given above for f € [fo, fi], a band of strictly 
positive frequencies, and is infinite for frequencies outside this band. 
Let b; = min; 6;. Then, for small crosstalk, 


Amin exp l(bi.f6 a cr fo), 
and 


logio Amin rz 0.434bi fi( 1 a Set ye 
k 
Thus, for small P, we have from (45), 


ee See 0.434, fi(1 — S28). (47) 
Al b; 


As in (46), the effect of the crosstalk is to reduce the dependence of 
channel capacity on cable length. 


4.3 Numerical example 


We consider a two-twisted-pair cable with white additive noise. The 
transfer function matrix H(f) is given by 


1 wrkli f 
=— —yl acre 
ce add ace 1 | 


with J in feet and f in hertz and 
y = av2rfi + ib2rf. 


The off-diagonal terms in the matrix H(f) represent far-end cross- 
talk. This model is an approximate representation of an experimental 
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Fig. 1—Channel capacity for experimental cable. Capacity C in units of 108 
bits per second is plotted as a function of signal-to-noise ratio P* for various values of 
cable length J. 


two-pair cable. Parameters obtained from measurement are 


k = 1.26 X 10-22, 
a = 0.23 X 107, 
6 = 1.48 X 107°. 
This model is valid in the range 10° S$ 1 S 50 X 10? feet, and 10°/2 
< f= 10’ Hz. 
Since the noise is assumed white, the ;’s are the eigenvalues of 
(H*H)— and are given by 
‘ ‘i a: exp(2V27 af?) 
M =p = At = EY 
1 + (27)? f?kl 
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The capacity equations (39) become 


f 
a | “sae Ore = Oa SP (48a) 
fo 
and 
fi K* 
C2 max { 0, loge ey df, (48b) 
fo 


where fp = 10°/2 and f; = 107. 

Numerical evaluation of (48) for various values of P and / has been 
performed and the results are given in Figs. 1 and 2 and Table I. The 
figures show C' vs P* for various values of 1. The C axis is linear in 


109 


108 


£ = CABLE LENGTH IN 
THOUSANDS OF FEET 


107 


10° 
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Vig. 2—Channel capacity for experimental cable. Capacity C in bits per second is 


plotted as a function of signal-to-noise ratio P* for various values of cable length l. 
The C axis is logarithmic. 
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Table | — Values of channel capacity C in bits per second for various values of signal-to-noise ratio P* 
and cable length. Exponential notation is employed for C (aEb =a x 10°). 








Z9L 4 AYOWAW HLIM TANNVHO NVISSAVS 


Signal- Cable Length / in Thousands of Feet 
to-Noise 
Ratio 
P* in dB 2 4 6 8 10 12 14 16 18 20 22 24 
—46 1.338E2 | 2.60E1 | 5.07E0 — —— — — - — — — — 
— 44. 2.10E2 | 4.12K1 | 8.180 | 1.57E0 — — — — — — — 
—42 3.32EK2 | 6.538E1 | 1.28E1 | 2.538E0 — — — — - _- — — 
—40 5.25K2 | 1.0812 | 2.08EK1 | 3.95E0 — — — — — - — — 
—38 8.29E2 | 1.68E2 | 3.21K1 | 6.29E0 | 1.21E0 — -—~- - —— — — — 
— 36 1.31E3 | 2.58E2 | 5.11E1 | 9.92K0 | 1.97E0 — — — — -— — —— 
— 34 2.06E3 | 4.07E2 | 8.06E1 | 1.58E1 | 3.04E0 — — — — — — a 
—32 3.23H3 | 6.48E2 | 1.27E2 | 2.51E1 | 4.94E0 — — — — — oo — 
—30 5.07E3 | 1.01E3 | 2.01E2 | 3.96E1 | 7.86E0 | 1.52E0 —- — — — ~~ — 
—28 7.9413 | 1.59E3 | 3.17E2 | 6.27E1 | 1.2461 | 2.36E0 —- — — — —- — 
— 26 1.24K4 | 2.50E38 | 4.99E2 | 9.92E1 | 1.95F1 ) 3.79E0 —~ — ~- — — -— 
—24 1.938E4 | 3.91E3 | 7.87E2 | 1.57E2 | 3.10E1 | 6.11E0 | 1.19E0 — — — — — 
—22 2.98EH4 | 6.09E3 | 1.2683 | 2.47E2 | 4.91E1 | 9.66E0 | 1.89F0 —- — — — — 
—20 4.5914 | 9.47H3 | 1.94E3 | 3.90E2 | 7.72E1 | 1.52E1 | 3.030 — — — — — 
—18 7.0214 | 1.4714 | 3.08E3 | 6.14E2 | 1.22E2 | 2.41E1 | 4.77E0 —- — —- — 
—16 1.07E5 | 2.264 | 4.72E3 | 9.68E2 | 1.92E2 | 3.82E1 | 7.57E0 | 1.47E0 — — — — 
—14 1.60E5 | 3.45K4 | 7.838E3 | 1.51E3 | 3.083E2 | 6.00E1 | 1.18E1 | 2.29E0 — — — — 
—12 2.39H5 | 5.24H4 | 1.18H4 | 2.36E3 | 4.78E2 | 9.50F1 | 1.87E1 | 3.62E0 -— — os — 
—10 3.038H5 | 7.88E4 | 1.74E4 | 3.683 | 7.50E2 | 1.51E2 | 2.98E1 | 5.87E0 | 1.110 — — — 
—§ 5.14K5 | 1L17E5 | 2.65E4 | 5.70E3 | 1.18E3 | 2.36E2 | 4.73E1 | 9.16E0 | 1.77E0 — — — 
—6 7.40E5 | 1.73E5 | 4.0184 | 8.80E3 | 1.843 | 3.72E2 | 7.42E1 | 1.48E1 | 2.88E0 _— — — 
—4 1.05EK6 | 2.51E5 | 6.01E4 | 1.8584 | 2.87E3 | 5.87E2 | 1.17E2 | 2.30F1 | 4.45E0 — — — 
—2 1.47E6 | 3.60E5 | 8.91E4 | 2.06E4 | 4.4513 | 9.18E2 | 1.84E2 | 3.67E1 | 7.32E0 | 1.350 — — 
0 2.03E6 | 5.10E5 | 1.31E5 | 3.11E4 | 6.87E3 | 1.44B3 | 2.91E2 | 5.97E1 | 1.15E1 | 2.25E0 — — 
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Table | — continued 
Signal- Cable Length / in Thousands of Feet 
to-Noise 
Ratio 
P* in dB 2 4 6 8 10 12 14 16 18 20 22 24 


| ES | ES | | | | NS | A | A | | | 


2 2.76E6 | 7.10E5 | 1.89E5 | 4.66E4 | 1.05E4 | 2.24E3 | 4.5862 | 9.12E1 | 1.81E1 | 3.56H0 = = 

4 3.71E6 | 9.74E5 | 2.69E5 | 6.90H4 | 1.61H4 | 3.48H38 | 7.192 | 1.44E2 | 2.88E1 | 5.610 |; 1.110 x 

6 4.90E6 | 1.32E6 | 3.78E5 | 1.015 | 2.48E4 | 5.388E3 | 1.12E3 | 2.27K2 | 4.55E1 | 8.99K0 | 1.74E0 oe 

8 6.40E6 | 1.75E6 |} 5.23E5 | 1.465 | 3.644 | 8.263 | 1.75E3 | 3.582 | 7.12E1 |} 1.42E1 | 2.800 = 
10 8.24E6 | 2.29E6 | 7.11E5 | 2.07E5 | 5.3894 | 1.26H4 | 2.7383 | 5.62E2 | 1.18E2 | 2.24K1 | 4.82K0 = 
12 1.05E7 | 2.966 | 9.538E5 | 2.90E5 | 7.894 | 1.914 | 4.21K3 | 8.79E2 | 1.78H2 | 3.541 | 6.95E0 | 1.360 
14 1.32E7 | 3.78E6 | 1.26E6 | 4.00E5 | 1.14E5 | 2.86E4 | 6.48E3 | 1.87E3 | 2.802 | 5.54E1 | 1.09E1 | 2.200 
16 1.6467 | 4.766 | 1.68E6 | 5.48E5 | 1.62E5 | 4.25E4 | 9.91E3 | 2.13E3 | 4.38E2 | 8.79E1 | 1.73E1 | 3.440 
18 2.02E7 | 5.92E6 | 2.09E6 | 7.26E5 | 2.27E5 | 6.23E4 | 1.50EK4 | 3.3183 | 6.88E2 | 1.39E2 | 2.73E1 | 5.430 
20 2.46E7 | 7.27E6 | 2.64E6 | 9.538E5 | 3.145 | 9.01E4 | 2.26H4 | 5.09E3 | 1.07E3 | 2.19E2 | 4.351 | 8.490 
22 2.98E7 | 8.85E6 | 3.29E6 | 1.23E6 | 4.26E5 | 1.285 | 3.36E4 | 7.80E3 | 1.67E3 | 3.42E2 | 6.83E1 | 1.37kK1 
24 3.57E7 | 1.07E7 | 4.06E6 | 1.58E6 | 5.685 | 1.805 | 4.944 | 1.19E4 | 2.60E3 | 5.388E2 | 1.082 | 2.12K1 
26 4.25E7 | 1.27E7 | 4.95E6 | 1.98E6 | 7.47E5 | 2.49E5 | 7.16E4 | 1.79E4 | 4.00E3 | 8.41E2 | 1.71E2 | 3.40E1 
28 5.02E7 | 1.5187 | 5.97E6 | 2.46E6 | 9.66E5 | 3.39E5 | 1.025 | 2.66H4 | 6.14E3 | 1.31E3 | 2.682 | 5.34E1 
30 5.8817 | 1.78E7 | 7.13E6 | 3.02E6 | 1.28E6 | 4.53E5 | 1.44E5 | 3.92E4 | 9.36E3 | 2.04H3 | 4.21E2 | 6.49K1 
32 6.86E7 | 2.07E7 | 8.44E6 | 3.67E6 | 1.556 | 5.96E5 | 2.00E5 | 5.714 | 1.41B4 | 3.153 | 6.59H2 | 1.3382 
34 7.92E7 | 2.40E7 | 9.92E6 | 4.41E6 | 1.92E6 | 7.72E5 | 2.72E5 | 8.19E4 | 2.11E4 | 4.8463 | 1.03E38 | 2.10H2 
36 9.04E7 | 2.77E7 | 1.16E7 | 5.25E6 | 2.36E6 | 9.86E5 | 3.65E5 | 1.16E5 | 3.13E4 | 7.39E3 | 1.60E3 | 3.29h2 
38 1.02E8 | 3.18E7 | 1.8467 | 6.19E6 | 2.86E6 | 1.246 | 4.82E5 | 1.61E5 | 4.57H4 | 1.12E4 | 2.48H3 | 5.1512 
40 1.14E8 | 3.62E7 | 1.52E7 | 7.25E6 | 3.43E6 | 1.54E6 | 6.26E5 | 2.21E5 | 6.584 | 1.68E4 | 3.82E3 | 8.06E2 
42 1.26E8 | 4.11E7 | 1.76E7 | 8.42E6 | 4.07E6 | 1.89E6 | 8.01E5 | 2.97E5 | 9.344 | 2.49E4 | 5.8463 | 1.26K3 
44 1.39E8 | 4.64E7 | 2.01E7 | 9.71E6 | 4.79E6 | 2.29E6 | 1.01E6 | 3.98E5 | 1.31E5 | 3.664 | 8.88E3 | 1.95E3 
46 1.51E8 | 5.22E7 | 2.27E7 | 1.11E7 | 5.60E6 | 2.75E6 | 1.266 | 5.12E5 | 1.79E5 | 5.29E4 | 1.33E4 | 3.01K3 
48 1.68E8 | 5.857 | 2.56E7 | 1.27E7 | 6.50E6 | 3.26E6 | 1.54E6 | 6.58E5 | 2.43E5 | 7.544 | 1.99E4 | 4.613 
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Table |— continued 


Signal- Cable Length / in Thousands of Feet 
to-Noise 
Ratio 
P*in dB 2 4 6 8 10 12 14 16 18 20 22 24 


—— | | |S | |S | fc | Se |S fn | See | sienna a 


Tig. 1, and is logarithmic in Fig. 2. The linear regions discussed above 
are evident in these figures. The asymptotic estimates for large P 
given in (41) and (46) are AC/AP* 6.3 X 10° (b/s/dB), and 
AC/Al = — 70 X 103 (b/s/ft). For constant C, AP*/Al% 11 x 1073 
(dB/ft) or about 58 dB/mile. This is the amount of increase of signal- 
to-noise ratio necessary to maintain a fixed level of C as length is 
increased. The asymptotic estimates for small P given (48) and (47) 
are A logyC/AP* 1/10, and A logyw:C/Al & — 0.353 X 107%. For 
constant C, AP*/Al& 3.53 X 1073(dB/ft), or about 19 dB/mile. 
These asymptotic estimates are borne out in the numerical evaluation. 
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APPENDIX 
A.1 Proof of Theorem 2 


a. Let y(n) = (Fx)(n) = Li f(n — k)x(k), and let 2d, || f(n) 
=C < o,. From the triangle and Schwarz inequalities, 


lye) I? SC Fe — Bll ey)? Ss CX MN fle — Ile). 


Hence, 
liv? = Eyl? SCX E Ie — Hl-leG|? Ss CAllall 
and |s| < C. 


b. From Parseval’s theorem, 


,  f IF@X@llas 
= 27 _______. = max ||F'(6)||?. 


{| |x|]? [ |X (8) ||2d0 8 


a 
ald 





ce. If F is self-adjoint, then 
(x, §x)| S [S| -[[]x[]P S max []'(6)| -IIIxlil?. 
d. The essence of this result is a matricized version of Wiener’s 
well-known theorem on the reciprocal of an absolutely convergent 


Iourier series (see Ref. 5, p. 430). 
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Suppose det (6) + 0. We show that § has a bounded inverse in L. 
Now det F(@) is a scalar function consisting of a sum of products of 
functions each possessing an absolutely convergent Fourier series." 
Hence, det F(6) and, by Wiener’s theorem, [det F(@) ]-! have abso- 
lutely convergent Fourier series. Each element of F-'(@) is the ratio 
of a minor determinant to det /’(6@) so each element has an absolutely 
convergent Fourier series. Hence, #’~1(@) hasa Fourier series >, n g(n)e*”® 
with >>, ||g(n)|| < ©. Consequently, 5 has a bounded inverse 37} 
in L. 

Conversely, let § have a bounded inverse $1. Then there exists 
an a> O such that for all x in I (— ~, ©), |{|x{|| = al||xi]]. 
Let X(@) = >) n x(n)e*”®. From Parseval’s theorem, 


—T 


[7 \FP@x@\Pa8 


0<a?< inf - 
“Of x@|Pa0 


which implies, since /(@) is continuous, that /(@) is one-to-one; i.e., 
det F(6) # Ofor —7 S 6S zz. This completes the proof of Theorem 2. 

There are other interesting properties of the class L, which are not 
directly relevant to the main results of this paper. For the sake of 
completeness, we mention two generalizations of Theorem 2d, that 
also have well-known scalar counterparts: 


(2) Let § € L and let o (F) denote the set of all eigenvalues of 
F(6), —x7 S607. Let I be the identity on Jf. For \ any 
complex number, AJ — §& has a bounded inverse in L if and 
only if \€ a(S). 

(iz) Let g(-) be any function analytic in a neighborhood containing 
a(S). Then there is an operator g(S) in LZ which has as its 
transfer matrix the function g[ F(@) ]. 


A.2 Lemmas 3 and 4 


These lemmas apply to the special case where (0) = I;, and 
R(6) = R2(@). Let x, y, and z be the channel input, output, and noise 
sequences respectively, and let x“, y‘%), z\) be the corresponding 
finite sequences. Thus, 


yO) = x) 4 20), 
* Note that Sallf(n) ||< © if and only if Sialfij(n) |< forl S$i,j Ss. 
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Letting T'y be the whitening matrix defined in (18), we have 
V 7 Try = T yx) -- W, 


where Eww! = Iy.,. This is precisely the discrete-time version of the 
problem treated in Chapter 8 of Gallager.! The results obtained there 
apply here exactly when we use, instead of his Lemma 8.5.2 (the 
Kac-Murdock-Szego theorem), the following discrete-time version: 


Theorem. Let {c;} 7 = 0, +1, --- be a sequence of s X s matrices such 
that the Ns X Ns matrix Cy = {c:-;}, 1, 7 = 0, ---, N — 1, ts Hermit- 
tan, and dix ||cx|] < ©. Let vf, vf, ---, o&P be the eigenvalues of Cy 
(each counted according to its multiplicity) and let d1(@), \2(8), ---, As(8) 
be the eigenvalues of C(@) = >> x cxe**®. Let g(-) be any continuous function 
defined on an interval containing the values {\,(0): —7 SOS 7,k = 1, 
2, +++, s}. Then 


: 1 sn ww) — 8 7 
tim 57 X90) = 5 Xf oLns(0) Jao, 


Furthermore, let Dy(x) = 1/N (number of eigenvalues vf < x). Then 


im Dos f dé. 
N—-« 2m Kat Jd4(0) Sz 
In the scalar case, s = 1, this theorem represents well-known results, § 
a simple account of which can also be found in Ref. 7. The validity of 
the theorem for s > 1 follows on verifying that the arguments em- 
ployed in Ref. 7 are valid in general for s = 1. 


A.3 Derivation of (38) 


We show how to obtain the capacity formula given in (38). We will 
do this for the scalar case; the result for vectors follows similarly. The 
capacity formula justified in Theorem 1 can be stated as follows. 

The channel input and output are sequences x = {z,}2. and 
y = {yn}2. respectively, related by 


Yn = 2 heb + ny (49) 


k=—o 


where h = {ha}lna=—. is a fixed sequence and z = {z,}*%. is astationary 
sequence of Gaussian random variables for which Ez, = 0, 


Hemem—n = Pay ~2< nm< @®, (50) 
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We write (49) symbolically as 
y = hax +2, (51) 


where ‘‘*”’ denotes vector convolution. The capacity of this channel is 
given as follows. 

For codes of block length N, say that the code vectors x must 
satisfy 


tr =0, nE€El[0, N — 1]. (52a) 
TS es 2b 
N 2a in = OD. (5 ) 
The capacity is then given by 
ae Or ave). ) 
C=r = df max (0, loge 7 ) (53a) 
where 
Ho(f) = XS hes, | f| < W, (53b) 
and 
Ro(f) = X Pret ns, | f| SW, (53¢) 


are the discrete Fourier transforms of {h,} and {7,} respectively, and 
kK is the unique solution of 


So = az [di max (0, K -— 2, ). (53d) 


A.3.1 Some facts about band-limited functions 


Before discussing the continuous-time channel, we digress to mention 
several facts about band-limited functions. 

We denote time functions by lower-case letters, e.g., u(é), and the 
corresponding Fourier transform by upper-case letters, e.g., U(f). 
Thus 


io a) 


U(f) = i, _ (ded, (54a) 


and 
2 [ U(fye-2"tdf. (54b) 


We shall assume that all functions are square-integrable, and all 
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integrals and infinite sums are limits in the mean. We say that w is 
band-limited to W Hz if U(f) = 0, |f| > W. Let 


. n 
sin 2rW (: — ai) 


gn(t) = eerey 


be the sampling functions. Note that for k,n = 0, £1, -:-, 


k\ |0 k#n, 
o( sy) = aw ean eee) 


n=0,+1,42,-:-, (55) 


and 
elinw/W)f 


: fl <W 
= i2af = 
Ga(f) [- gn(te2*Ftdt | : > Ww (56b) 
(so that g, is band-limited), and for k,n = 0, +1, --- 
[aaa = [" acnexnas 
0, n#Fk, 


W 
a (irf/W)(n—k) — 
Je af a ees (56c) 


Further, the well-known Sampling Theorem implies that any band- 
limited function, w(t), can be written as 


ut) = DY ungn(d) (67a) 


where 


] n 


Further, from (57), and (56c), 
J ” wdt=2W Xt ue. (7c) 


Let u(é) = © vagn(t) and v(t) = > vagn(t) be band-limited functions. 
Then their convolution is | 


wi) = fo ult — oat = X wagald), (58a) 
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where 


ea] 


i= DD Mela SOO SS Os, (58b) 


1.e., W = sv. 
Finally, let z(t) (— ~<t< «) bearandom process with Ez(t) = 0, 
and covariance 


EeQeim~ r) = rir) — eo <tir< o. 
Let 


R(f) = [ r(t)ei* dt 


satisfy R(f) = 0, |f| > W, so that r(t) is band-limited. Then, from 


(57), 
HO = Lary ( gap ) a0 
and from (56b), 


ROD = Cage (gap) GoD = gp De ( gap Jeers. 60) 


Further, the random process z(t) is a band-limited function so that by 
(57) we can write 


al) = L zagald) = X gap? ( gop ) on: 
Thus, 


fa © Henan ay ®( air) ("a") 7 aire" (air) 


Thus, the discrete Fourier transform of {7,,} is, using — 
: 1 
= p ol(ir/W)nf — — (ir/W)nf — 
Ro(f) os P ne aw): pa (sir) e rR f). (60) 


A.3.2 The continuous-time channel 


The continuous-time channel is defined as follows. The channel in- 
put and output are functions x(t) and y(¢), respectively, where 


y(t) = i ht — r)a(\) dd + a(t), -2 <t<, (61) 
where /A(é) is a fixed function and 2(¢) is a Gaussian random process 
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with covariance as described above. We assume that 2x(¢), h(t), and, 
therefore, y(t) are band-limited to W Hz. Let us expand a, h, z, and y 
into series in g,(t) as in (57). Using (58) we obtain 


Yn = Yo ha—mBm + 2m- (62) 
Since knowledge of the sequences of coefficients {z,}, {yn}, ete. is 
equivalent to knowledge of the time functions xz(f), y(t), etc., the 
continuous-time channel is equivalent to the discrete-time channel 
discussed at the beginning of this appendix. It remains to find the 
corresponding parameters. 

Now the code words (in a J-second block-coding interval) are 
taken to be band-limited functions z(t) such that the samples 
a(n/2W) = 0, for (n/2W) <0 or (n/2W) 2 T, i.e., x(n/2W) = 0, 
n€&[0,N — 1], where N = 2WT. Thus, x, = (1/2W)z(n/2W) = 0, 
n€[0, N — 1]. The condition 


J ” 2%(t)dt < ST 
is, in the light of (57c), 


Let 2 s 
No QW)? 


IIA 


(63) 





The quantity 
Ho(f) = dX bret! = ThGalf) = H(f), 


and by (60), Ro(f) = (1/2W)R(f). Thus, the continuous-time channel 

is equivalent to a discrete-time channel with Sp = S/2W, Hp(f) 

= H(f), and Rp(f) = (1/2W)R(f). Thus, from (53) 

ee 
fi(f) 


) 


(a2 57 a 
- ay |. f max (0, O22 


where K is the solution to 


S ee W oe 0 R(f) 
ai = aw |, oem (0, * — 37 TH) 3) 


Letting K* = 2WK, we have 


S= f afmax(0, K* - HOR 5) 
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and 


C= ap [af max (0, log, AWE), 


The expression for C is in bits per sampling time. To obtain C in 
bits per second, multiply by 2W. 


A.4 Proof of equation (42) 


Equation (42) follows at once from the theorem below if the func- 
tions \j,7 = 1, 2, ---s, are replaced by a single function f representing 
a concatenation of the functions \;; f will be bounded away from zero 
since the same is true of each )j. 


Theorem: Let f(-) be a measurable function on a finite interval and 
ess inf f(x) = fo>0. For any K > fo define A = {2; f(x) S K} 
and let 6 be the measure of A. Define 


log 5 [ Ot eee ah log f(a) dx 


Pi = 
Kise sf f(a)dx 


Then 
lim [,(k) = 0. 


K->fo 


Proof: Without loss of generality we can take the log to be the natural 
log and can assume fo = 1. Let f= 1-+gand K =1+k. For each 


n = 1 define 
g? = ;/ g™(x)dx. 
OJaA 


For all cx € A, g(x) Sk and g Sk. For k < 1, the log may be ex- 
panded in a power series. After some rearrangement of terms, we have 


fo) 


L a a 

os, Qn ~— Fan\) — Qn+1 _. A2n+1 
5 la @ ~~ a | 
1,(K) = = . 





But g S k and by Jensen’s inequality g” S g” for all n = 1. Then 


2n — Fen 
nt g’") 
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IV 
= 


Now for all n 


n—1 
kn — gr = (bg) & kehig! < (bk — g)nkr, 
1=0 


and 


k 


Krol 12 
22 = 
2» I ] ae k2? 


(XK) S$ SO +9") s 





a] 


or, since K = 1+ &, 


kK — 1 
[ (Kk) s > 
r( STH (ke 1)?’ 
and the result is proved. 
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Experimental Results on a Single-Material 
Optical Fiber 


By R. D. STANDLEY and W. S. HOLDEN 
(Manuscript received November 20, 1973) 


Experiments confirm a number of theoretical predictions regarding the 
behavior of single-material optical fibers. In particular, the predicted modal 
velocity spread (10's of ns/km) and the numerical aperture are found 
to be 1n good agreement with theory. 


I. INTRODUCTION 


The concept of optical fibers made from a single material was the 
subject of a recent article;! the advantages of such fibers are their 
construction of low-loss fused silica and their freedom from the 
problems associated with glass interfaces. In this paper, we present 
experimental results on the dispersion observed in a particular single- 
material fiber. The significance of this work lies in the good agreement 
between the theory and the experiment; the results do not represent a 
careful evaluation of single-material fiber as a transmission medium. 

The results of the experiments are as follows: 


(2) The predicted modal velocity spread was confirmed. 
(22) The measured numerical aperture is directly proportional to 
wavelength, as predicted by theory. 
(4272) There was very little mode coupling between lowest order 
modes for lengths up to 100 meters. 
(iv) Penetration of energy into the support structure was small (the 
decay constant is about 20 dB/yum). 


Il. THEORY OF MODAL VELOCITY SPREAD 


The theory of single-material fibers is summarized in Ref. 1. In 
the present discussion, we concentrate on the modal dispersion of such 
fibers. Given the structure shown in Fig. 1, we assumed that the 
electromagnetic fields vary either sinusoidally or exponentially along 

779 





a 


Fig. 1—Idealized fiber cross section. 


x, y, and z. We further assume the modes of interest to be far above 
cutoff so that the transverse wave numbers can be considered to be 
constant; we approximate them as 


ets 

k= A (1) 
and 

hy = (2) 


A and B are defined in Fig. 1 and yp and » are the mode orders. It then 
can be shown that the group velocity of the mode of order uv is ap- 


proximately 
2 2 2 
Vow =$ 1-3 [(4) +(#) |} 8) 


This gives for the expected time spread between the fastest and the up 


pulses 
_ LX uw \2 y \2 
sro (G) + (a) | ' 


where Z is the fiber length. From Ref. 1, the maximum value of the 
bracketed term in eq. 4 is b-; b is defined in Fig. 1. Thus, the max- 
imum time spread between pulses is 


AtTmax = 2 Gy: (5) 


Sen 
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Ii. EXPERIMENTAL RESULTS 
3.1 Description of the fibers 


The fibers used for our experiments were made of fused silica and 
had the cross section shown in Fig. 2. Because of the multimode 
nature of this particular fiber, we anticipate that the theory developed 
in Section II for the rectangular fiber geometry will still approximately 
apply. 

The following numerical result can serve as a guide in predicting 
the actual measurement. Assume a fiber with A = B= 10 um, 
b = 2 um, and n = 1.46, and let \ = 1 wm. Then, 


Atay = 2.85L(u? + v?)ps. (6) 


For a 100-meter fiber length, the low-order pulses would have time 
separations of a few tenths of a nanosecond. This implies the following: 
if such a fiber were excited by a pulse whose width is small compared 
to Ar,z,, then, assuming little mode coupling, each mode should be 





Fig. 2—Photograph of experimental fiber cross section. 
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Fig. 3—Fiber output as a function of time and launching conditions. 
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output as several pulses each representing a mode group. (In this case 
the fiber was about 90 meters long.) The relative amplitudes of the 
output pulses could be changed by varying the input launching 
condition, which is demonstrated by the two photographs shown in 
Fig. 3. 

Additional information was obtained by measuring the velocity 
difference between the output pulses as a function of fiber length; 
Fig. 4 shows the results for the first few modes. The results are in 
reasonable agreement with that predicted by eq. 4; the predicted time 
spread between the first four modes is 8 ps/m, 17 ps/m, and 21 ps/m, 
whereas the measured values were 7 ps/m, 16 ps/m, and 22 ps/m. 
Note that the time difference approaches zero for zero length, which 
suggests little mode coupling among the lower-order modes as does the 
previous observation that individual modes could be preferentially 
excited. Although not shown, the higher-order modes behave differ- 
ently, having a time difference which appears to be related to the lower 
order modes. This suggests that the higher-order modes are possibly 
generated by mode coupling from the lower-order modes and that the 
higher-order modes suffer more loss. 


RELATIVE MODAL DELAY IN NANOSECONDS 





0 10 20 30 40 50 60 70 80 90 100 
FIBER LENGTH IN METERS 


Fig. 4—Relative modal delay as a function of fiber length. 
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3.3 Numerical aperture 


Theory predicts a numerical aperture (half-angle of the fiber radia- 
tion cone) of NA = \/2b. For our fiber, b was 2 um so that the theoret- 
ical NA at\ = 1.06 um was 0.265, and at \ = 0.6328 um the NA should 
be 0.1582. The predictions compare favorably with our measured re- 
sults, which were 0.27 at \ = 1.06 um and 0.16 at \ = 0.6328 um as 
determined from the angular spread of the far field radiation. 


IV. SUMMARY 


The experimental investigation of the dispersion in a single-material 
optical fiber showed good agreement with theory. The linear depen- 
dence of numerical aperture on wavelength also was verified. 
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Optimal Trade-Off of Mode-Mixing Optical 
Filtering and Index Difference in Digital 
Fiber Optic Communication Systems 


By S. D. PERSONICK 


(Manuscript received December 11, 1973) 


In a digital fiber optical communication system, the optical power re- 
quired at the receiver input to achieve a desired error rate depends upon 
the shape of the received pulses. In systems employing multimode fibers 
and/or broadband sources, we can experience pulse spreading in propaga- 
tion because of the group velocity differences of different modes or because 
of dispersion. In an effort to control or compensate for pulse spreading, 
we can trade off coupling efficiency between the light source and the fiber 
(by varying the core-cladding index difference or bandlimiting the source), 
scattering loss in the fiber (by introducing mode coupling), and equalization 
in the receiver at baseband. This paper investigates the optimal trade-off 
for various fiber-source combinations. 


I. INTRODUCTION AND REVIEW OF BACKGROUND MATERIAL 


In digital fiber optic communication systems, as in other digital 
systems, the received power required at a repeater to achieve a de- 
sired error rate depends upon the shape of the received pulses. A 
previous paper! showed that the minimum average power requirement 
results from a pulse that is sufficiently narrow so that its energy spec- 
trum is almost constant for all frequencies passed by the receiver 
(ideally, an impulse). For other received pulse shapes, we can define 
the additional power required, in decibels, as a ‘power penalty” for 
not having impulse-shaped pulses. Typical calculations of this power 
penalty for ‘on-off’ signaling and a receiver employing avalanche 
gain with a high impedance front end! are shown for various families 
of received pulse shapes in Fig. 1. In that figure, the parameter o/T 
is defined as follows: 


a= Fi [ro@vat — ral hetteat (1) 
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EXPONENTIAL 
(AL) 


POWER PENALTY IN DECIBELS 


~____ RECTANGLE 
(LU) 





o/T 


Fig. 1—Typical power penalty vs o/T. 


where 
T = spacing in time between binary digits 

hp(t) 

A 


I 


received optical pulse shape 


area under /A>(t). 


We shall refer to o as the rms pulse width.* 

It has been shown? that, in long fibers, the ‘power impulse re- 
sponse”’ of the fiber approaches a Gaussian shape. In the rest of this 
paper, we assume that the received optical pulse is Gaussian in shape 
and that it has an rms width determined by the fiber delay distortion. 
That is, we assume that the rms width of the fiber input pulse is 
sufficiently small so that the power penalty associated with the rms 


“In a Gaussian-shaped pulse, the rms width is about 0.425 times the full width 


ee Aa ea points. In a rectangular pulse, the rms width is 1/12 the 
ull width. 
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sum of that width and the rms fiber impulse response width is the same 
as the power penalty associated with the rms fiber impulse response 
width alone. 

From various heuristic analyses (see the appendix), we can conclude 
that the rms width of the received pulse is approximately the rms sum 
of the delay distortion in the fiber resulting from material dispersion 
(because of the variation of group velocity with wavelength associated 
with the use of a broadband source) and the delay distortion associated 
with the spread in the group delays of various fiber modes (when a 
multimode fiber is used). That 1s, 


G = (Oaiinerdon ape. nda) (2) 


where Guispersion ~ Optical bandwidth: fiber length = B-Z and where 
Omode 18 determined as follows’ 


Case 1. Conventional clad multimode fibers without mode coupling. 


Fin = 0080 = Ee 
C 


where 

n = index of refraction of the core 

A = (index of refraction of the core — index of refraction of the 
cladding) /n 

c = speed of light. 


Case 2. Conventional clad fibers with complete mode coupling after 
a distance Lc. 


Smoae = 0.289 on Vibe toe Este. 

= 0.289 1 for L < Le. 

Case 3. Ideal graded index fiber without mode coupling. 
Coen = 0.087 ** Fr. 


Case 4. Other fibers? can be treated once the techniques outlined below 
are understood. 


From the definitions of caispersion 20d Cmoae ADOVe, we see that, for a 
broadband (incoherent) source, the material dispersion contribution 
to « can be controlled by limiting the optical bandwidth B being used. 
However, if the optical source bandwidth, B., must be reduced by 
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filtering, then the average power into the guide will be reduced by the 
factor B/B,. Similarly, we can control the mode delay spread by reduc- 
ing the index difference A. If a multimode (incoherent) source is being 
used, the average power into the guide is proportional to A. Thus, we 
trade off input power against mode delay spread. Furthermore, even 
when we use a coherent source which in principle can be focused into 
any fiber, we must be careful when A becomes significantly smaller than 
0.005, since the fiber loss at bends becomes large. (Exactly what value 
of A is too small to be practical 1s an open question.) For fibers with 
mode coupling, we can control omode by decreasing Le (increasing the 
mode coupling). However, this causes the coupling to radiating modes 
to increase, thus increasing the fiber loss.’ Here again, there is a trade- 
off between the average power we receive at the fiber output and the 
received rms pulse width. For a fixed shape of the mechanical spectrum 
of fiber geometry perturbations, the radiation loss per unit length of 
the fiber resulting from mode coupling is inversely proportional to Le, 
1.¢@., 

radiation loss in nepers = a,L/Le, (3) 
where 


a = constant depending upon the shape of the mechanical spectrum 
of the geometry perturbations causing coupling (and possibly 
upon the index difference A). Lc depends upon the amplitude of 
the mechanical perturbations. 


In the following sections we derive the optimal trade-off between a, 
B, A, and Le for various combinations of sources and fibers to maximize 
the allowable fiber length ZL between the optical source and the 
repeater. 


ll. ANALYSIS 
2.1 Incoherent source, conventional clad fiber, no mode coupling 


Let the average power into the guide be P, when the index difference 
A is at some maximum practical value A, and when the full source 
optical bandwidth B, is being used. Let the loss of the fiber be a nepers 
per kilometer. Let the power penalty from the nonzero value of the 
received rms pulse width, in nepers, be f(¢/T). Let the required power 
at the receiver be P, when ¢/T = 0. If we use a value of A S A, and 
filter the source output to have an optical bandwidth B < B,, then we 
must have. 


= 


P, ; eth > PreftolT) (4) 


Plp 
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where, from (2), 


2 2) 4 
+= (Ba) + (2a) = mera 


and C, C2 are constants. 


We rewrite (4) as follows (using equality to maximize L), 


P AB a 
—S gab — |" — e-fiT)\ , 
eae 


r 


To maximize L, we must choose A/A, and B/B, to maximize the term 
in braces subject to the constraint that these ratios cannot exceed 
unity. We define —10 log (term in braces) as the ‘“‘excess loss.”’ 

By equating appropriate partial derivatives of the excess loss to 
zero, we obtain the following equations for optimizing A and B (A and B 
which minimize excess loss). 


2 
ig acts Pape A 
f ( 7) oe J provided A/A, S 1, (5a) 
otherwise A/A, = 1, 
2 
/ us Ga — 1 R 
ci ( > | i 1 provided B/B, < 1, (5b) 
otherwise B/B, = 1 


; 


where f’(z) = d/dz[ f(x) ]|,-, and om and og are defined in (4). For 
sufficiently long lengths L, where both A/A, and &/B, are less than 
unity, we obtain [by adding (5a) to 5(b) ] 


on = 04 = oe, 6) 


where « is the solution f’(x)z/2 = 1. More specifically, we obtain the 
following 
BL _ «T 


B, v2’ 


A. 








(7) 


S|5 


Cy C 
and therefore’ 
l l xT” —f (2) 
excess loss = —10 log seas e | 
for 


L > maximum of Be, a 
fe V2C1’ v2C2 


* Throughout this paper we use the parameter L/T (guide length/time slot width) 
frequently. The larger the fiber length or the smaller the time slot width, the more 
excess loss must be incurred to control or compensate for pulse spreading. 


DIGITAL FIBER OPTIC COMMUNICATION 789 


From the Gaussian power penalty curve of Fig. 1, we obtain the value 
of x where f’(x)x = 2 to be 0.37. At that value of z, 


—10 log e4%@ = 3.3 dB. 


As L/T decreases, either A/A, or B/B, will eventually reach unity. 
When that happens, o will approach og or om, respectively, for shorter 
lengths LZ. Then, from (5), ¢/T will approach the solution of f’(z)z = 1. 
Furthermore, the excess loss will approach either 


excess loss > —10 log Ee crn | (8) 
1 


if A/A, reaches unity first and L/T > z/C, or 
2T 
_ eo oy le) 
excess loss — — 10 log Ie, € (9) 
if B/B, reaches unity first and L/T > 2/Co. 


From the Gaussian curve of Fig. 1, the solution of f’(z) z= 1 is 
z = 0.3. At that value of z, —10 log e-/%) = 1.8 dB. 


7 
/ 


x 
START EQUALIZING / 


a“ 
START CUTTING 
INDEX DIFFERENCE 


FILTERING 


EXCESS LOSS IN DECIBELS 





10-3 10-2 10-1 109 
A 5 


Fig. 2—Iixcess loss vs L/T for conventional fiber. 
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For L/T sufficiently small so that both A/A, and B/B, = 1, the 
excess loss will approach zero as L/T decreases. 


Example 1 


In a fused silica fiber, the material dispersion has been measured at 
9 ps/km per angstrom of wavelength difference. With a typical GaAs 
LED* (light-emitting diode), this results in a value of og of 1.5 ns/km 
of pulse spreading at the full bandwidth B,. Thus, Ci can be set at 1.5 
ns/km. For a fiber with a maximum index difference A of 0.01, C; 
would be given by 14.5 ns/km. 

From (7) we obtain 

(0.377)? 

= BOE eg 5) (14.5) 


ef (0.37) 


excess loss 


for L/T > 0.174, where the length L is in kilometers and the time slot 
width 7' is in nanoseconds. 

From (8) we find that, for 0.0206 < L/T < 0.174, the excess loss 
asymptotically approaches 


2T T 
is Spotless ae ea 
excess loss — —10 log Fe € | 10 log T, + 18.6 dB 
for 0.0206 < L/T < 0.174. 

For L/T < 0.0206, the excess loss asymptotically approaches zero. 
Figure 2 is a plot of excess loss vs L/T in this example. 


2.2 Incoherent source, ideal graded index fiber, no mode coupling 


Following the same procedures as in 2.1, we can replace o» for a 
self-focusing fiber by C3(A/A,)?L (where A, is the maximum allowable 
value of A). We then obtain the following set of equations which 
determine the values of A and B that minimize the excess loss 


»{o\2on _ P 
7 GE: | provided om < 1, (Qa) 
otherwise A/A, = 1, 
2 
(eee ok oe B 
: (>) oT i provided Bi Be = 1, (9b) 
otherwise B/B, = 1. 


* Assuming a Gaussian-shaped optical spectrum with bandwidth between the half- 
power points of about 400 A. 
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For sufficiently long lengths L, where both A/A, and B/B, are less than 
unity, we obtain [by adding (9a) to twice (9b) | 


~ 


A \? x'T 
m= C:(=)L= 25 (10) 
B 2 
oq = CR L = x'T 3 
(x')#T? 4 \' 
excess loss = —10 log ee. eS (2") 
LIC WC; 


where x’ is the solution of f’(a’)z’ = 1.5, and where we must have 


L ; x’ 9 x fi 
— > maximum of } — . ]/% ey aes ee 
P | C1 V5 -" Cs V3 


For the Gaussian power penalty curve shown in Fig. 1, we have 
xz’ = 0.84 and —10 log e 4@ = 2.6 dB. 

As before, as L/T decreases, either B/B, or A/A, will reach unity. 
Thereafter, for smaller values of L/T we have the following: either 


aD 


g — oa; excess loss ~ —10 log Tc 
1 


ef), (11) 


vi, 
START / 
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START 
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Fig. 3—Excess loss vs L/T for graded index fiber. 
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where f’(z)z = 1, provided A/A, reaches unity first and L/T > z/C,, or 


g — om; excess loss > —10 log ( as yer | 
, LC; bf 

where f’(z2’)2z’ = 1, provided B/B, reaches unity first and L/T > 2’/Cs. 
For the Gaussian power penalty curve of Fig. 1, we obtain z = 0.3, 
—10 log e-%) = 1.8 dB, 2’ = 0.24, —10 log e-%@” = 0.98 dB. 

When L/T is sufficiently small so that A/A, and B/B, are both equal 
to unity, then the excess loss asymptotically approaches zero for 
smaller L/T. 


Example 2 


From Example 1 we have C, typically 1.5 ns/km. From (2) we have 
C; typically 0.019 ns/km for an ideal graded index fiber with a maxi- 
mum A = A, of 0.01. 

From (10) we obtain 


excess loss = —15 log ( z) + 4.85 dB 


for L/T > 10.8, where L is in kilometers and T is in nanoseconds. 
From (11) we obtain 


excess loss ~ — 10 log ( z) + 8.78 dB 


for 0.2 < L/T < 10.3. 

For L/T < 0.2, the excess loss asymptotically approaches zero as 
L/T approaches zero. | 

Figure 3 shows a plot of excess loss vs L/T for this example. 


2.3 Incoherent source, conventional fiber with adjustable mode coupling 


Now consider a conventional fiber with adjustable mode coupling. 
Using the same notation as in 2.1, we have the following condition 
from (2), (8), and (4)* 


Py ] i ls a i Mea LE (12) 


o Oo 


*As mentioned in Section I, the parameter a, depends upon the shape of the 
mechanical coupling spectrum and possibly upon the index difference A. Since this 
dependence upon A is not known analytically, we assume a, = constant independent 
of A. One could also assume a, « A% for some N (probably negative) and still obtain 
simple results similar to those that follow using analogous techniques. 
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where 
2 4 
p= (Cg 2) r (o.5 \EEe) } = a 


provided L > Lc. 
By setting appropriate partial derivatives to zero, we can minimize 
the excess loss given by 
excess loss = —10 log a= e—eoLhiLce—f(olT) \ . (13) 
Ay B, 
The optimizing values of A, B, and L¢ satisfy the following equations 


(we are assuming a, fixed by the shape of mechanical coupling 
spectrum). 








hid (5 om _ 1 provided A/A, < 1 
T }o es (14a) 
otherwise A/A, = 1, 
2 
a (ei ed es B 
f ( 4 oT 1, provided 6/B, < 1, (14b) 
otherwise B/B, = 1, 
al 4f¢ a2, 
i = f (7) aeT'? provided L> Le, (14¢c) 
= “0%, proydedaaetendt Sits. ida) 
C 


For sufficiently long fibers and if a < 0.5, we will have L > Le, 
A/Ao < 1, B/B, < 1, and therefore the following will hold :* 





excess loss = —10 log at i ae 1] (15) 
2 V2L?CiC Va. 
where x is the solution of f’(x)z/2 = 1 
aol 05 A _ eT B _ aT 
Le 7 Ao 20.DVa,’ = Bo SS VAL 1’ 
provided 
re odie ae rate ao S 0.5. 


It is convenient to consider L/T and L/Le as separate parameters. 


. *It is interesting to note that, with optimal mode coupling, in the region where 
A < Aj, the optimal value of A is increased by the factor V0.5/a, relative to the 
no-mode coupling case [see formulas for A/A, in (7) and (15) and also (8) and (16) ]. 


Further in this region (A < A,), the excess radiation loss from mode coupling (aol./L-) 
is always 0.5 neper. 
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As L/T decreases, either A/A, or B/B, will reach unity. 
If B/B, reaches unity first, then the excess loss will asymptotically 
approach the following for smaller values of L/T. 





2T 
a — =p Sf (2) 90.5 
excess loss > —10 log Cie e e (16) 
aly, A__ &f 
Le -_ Ao CoLy2q,’ 


where f’(z)z = 1, provided B/B, reaches unity first and 
L/T > 2/(C2 V 2a). 


If A/A, reaches unity first, then the excess loss will asymptotically 
approach the following for smaller values of L/T 


2T 
oe ——_ ef lz) 
excess loss — —10 log | iC, é (17) 
Be 
B, LC,’ 


provided A/A, reaches unity first and L/T 2 z/C1. 
For values L/T below that at which A/A, and B/B, both equal unity, 
the excess loss asymptotically approaches zero. 


Example 3 


Using the same parameter values as in Example 1 and assuming’ 
a, = 0.1, we obtain the following 


excess loss = —20 log - + 27 dB 


for L/T > 0.174 km/ns, 
excess loss => —10 log = + 17.3 dB 


for 0.0462 < L/T < 0.174, and 


excess loss = 0 for 7 < 0.0462. 


Figure 4 shows a plot of excess loss vs L/T for this example. 


* At this point, the achievable value of a, in practical fibers is a subject of specula- 
tion. We choose a, = 0.1 arbitrarily. 
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Fig. 4—Excess loss vs L/T for conventional fiber with coupling (a, = 0.1). 


2.4 Laser source, conventional fiber with adjustable mode coupling 


Here we assume that, to avoid excessive loss at bends, the index 
difference in the fiber, A, is fixed at some minimum allowable value 
Amm- To control pulse spreading we can trade off mode-coupling radia- 
tion loss against equalization penalty. The condition we must satisfy is 


Pe~ th e—eobllc = P,ef at), 


We wish to choose L¢ to minimize the excess loss given by 


excess loss = —10 log {e~eL/Lcg—-f\si/T)} (18) 
where* 
¢= CuvL Le, C, = 0.289 A nint/C. (19) 
We obtain the optimizing equation: 
—lisfo\a 
aoL/Le = 5 f (+)% (20) 


* Material dispersion is assumed negligible for a coherent source. That is, we assume 
single mode and a short-term bandwidth less than 1 A. 
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Fig. 5—Excess loss vs L/T—laser source, conventional fiber Amin = 0.001, a, = 0.1. 


To solve (20) we can pick a value of ¢/T and solve for f’(7/T) and 
f(c/T) graphically from Fig. 1. We then use those results in (20) 
to solve for L/Lc. Then we substitute into (19) to find L/T and into 
(18) to find the total excess loss. 


Example 4. 
Using Ay, = 0.001 and a, = 0.1, Fig. 5 shows a plot of excess loss 
vs L/T for this example. 


lil. APPLICATIONS 


If the optical power required at the receiver when the received pulses 
are very narrow is P, and the transmitted power (at maximum band- 
width and index difference) is P, and if the fiber loss in the absence of 
mode coupling loss is al, then we must have 


10 log (P,e-*”) — excess loss (ZL) = 10 log (P,) 


or, equivalently, * 


10 log < e--l = “excess gain (L)”’ = excess loss( LZ). 


* We define excess gain as number of decibels by which P,e~ *” exceeds P,. In effect, 
it is equal to the ‘‘allowable excess loss.”’ 
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At a given bit rate, 1/7, and given a, P,, and P,, we can plot excess 
loss (L) and ‘“‘excess gain (L)’’ simultaneously. The intersection of the 
two curves gives the maximum allowable distance LZ between the 
transmitter and the receiver. 


Example 5 


Assume that, at a bit rate of 25 Mb/s (7 = 40 ns), the required 
received power P, is approximately —58 dBm. Assume that a conven- 
tional fiber with mode coupling (a, = 0.1) and loss a = 5 dB/km is 
used. Assume that an incoherent source is being used with P, = —13 
dBm for A,x = 0.01. Assume that C; and C2 are 1.5 and 14.5 ns/km 
so that Fig. 4 applies. Figure 6 shows a plot of excess loss and excess 
gain vs L. It is apparent that the maximum length L between the trans- 
mitter and the receiver is 6.8 km. At that distance, the excess loss 
= 11.5 dB. Further, we have A = 0.0024, B/B, = 1 and Le = 1.36 
km.* 


IV. COMMENTS AND CONCLUSIONS 


The purpose of this paper has been to show how we can combine 
analytical results on fibers and repeaters to determine maximum re- 
peater spacings by optimization of available parameters. Since the 
fiber art is still young, many assumptions above are subject to ques- 
tion. We can summarize a few possible criticisms here. 

It is not known whether the assumption of the Gaussian pulse shape 
leads to overly conservative estimates of the equalization penalty. 
With time and experiments, the Gaussian pulse shape approximation 
will probably be improved upon. 

It is not known yet how much control the designer will have over 
the mode coupling and the index difference. Future analyses will 
have to take into account the practical constraints on these parameters, 

It is not known yet whether optical filters of the type assumed above 
can be built. Further, the above analysis neglects in-band insertion 
loss. 

It is hoped that, although the above analysis is somewhat simplistic, 
it can serve as a guide to the fiber system designer by pointing out the 
concepts and trade-offs involved. 


*It is interesting to note that, if mode coupling were not allowed, the excess loss 
curve would be about the same (calculated from Fig. 2) and therefore the maximum 


length L would be the same. However, at the maximum length A would be 0.0011. 
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Fig. 6—Excess loss and gain vs L. 


APPENDIX 


We wish to show heuristically that the total rms width of the fiber 
impulse response is the rms sum of the contribution from dispersion 
and the contribution resulting from mode delay spread. 

Suppose that, if the fiber is excited by a narrow-band source at wave- 
length d, the resultant output response is h(t). Let the mean arrival 
time and rms width of h(t) be defined as 


n= i / th (t) dt (21) 


1 


1 | 2 
7. = | a feroa — n| 


>» 


where 


Ay = [rat = area of h(t). 


Now suppose the fiber is excited by a narrow pulse from a broadband 
source having its output distributed in wavelength according to the 
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spectrum S(\). Intuitively,* we can write the fiber output response 
as follows: 


h(t) = | S(r)ha(t)dd = | rS(a) Aa] (22 ) dy. (22) 


Define $(A) as S(A)A, and let 


5 = [S@)Ardr = frou. 


Let S(A) = [S(\)A)/5 1]. It follows from (22) that the rms width of 
h(t) is given by 


| § [en 2 (5 foo) 


HT feSonaa] 
+| frasora — (f raySonar) |} (23) 


The first term in square brackets in (23) is a weighted average of the 
mean square width of the narrow-band pulse at different wavelengths. 
The second term is the mean square deviation of the narrow-band- 
mean-arrival time, i.e., the dispersion o%. If we next assume that 
o, & constant = o, (ie., that the rms width of the narrow-band im- 
pulse response is not dependent upon wavelength within the band of 
interest), then we obtain 


Q 
| 


g = lon + oj}, 


which is the desired result. 
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A Map Technique for Identifying 
Variables of Symmetry 


By L. M. GOODRICH 
(Manuscript received July 2, 1973) 


This paper presents a new map technique for identifying symmetrizable 
functions. The technique greatly reduces the work in ascertaining sym- 
metricity, and it 1s unique in being also applicable to completely or 
encompletely specified functions which: 


(2) Contain tmbedded symmetrizable function(s). 
(11) Are the complement of a function of type (1). 
(411) Contain an imbedded function of type (12). 


Discussion of the technique and its extensions is included. 


I. INTRODUCTION 


Recognition of symmetry in circuit design often can drastically 
reduce the problem of finding the least expensive circuit configuration. 
Multi-output circuits frequently have a symmetric circuit as a com- 
mon portion so no single error will result in a wrong output. 

As a consequence, numerous papers and chapters of books have 
presented recognition of symmetry in a switching function.1~* How- 
ever, none of these articles has presented a technique that is simple 
to apply and has natural extensions to accommodate both completely 
and incompletely specified functions that are almost symmetrizable.* 
This paper presents such a technique. 

Caldwell!? has demonstrated a technique using Karnaugh maps for 
recognizing symmetrizable functions (SF’s) of three or four variables, 
and has also demonstrated a procedure for extending this to functions 
of more variables. The extension requires the use of a large number of 
maps and the use of an expansion theorem a multiplicity of times. The 
Caldwell technique requires mapping all possible submaps in four of 
the variables. 


*See Section II for definitions. 
801 


Few authors have dealt with any subset of almost SF’s (ASF’s), 
and even fewer have tried to give practical examples of a technique 
for their solution. 

Born and Scidmore® have dealt with the special subset of ASF’s 
commonly called partial symmetric functions. They have solved a 
function in six variables for which, by their definition, a minimum 
solution exists and a resulting realization is shown in Fig. 1.* This same 
example has been solved by the technique presented in this paper and 
the resulting circuit, shown in Fig. 2, is significantly more economical.t 

In general, other techniques attempt to ascertain symmetry and find 
the center of symmetry (COS) simultaneously. The technique pre- 
sented here first uses a set of overlapping maps to ascertain what the 
COS must be if the function is an SF (or ASF) and then verifies whether 
the function is symmetric (or almost symmetric) about that COS. 


ll. METHOD 
2.1 Theory 


Shannon first stated the definition of a symmetric function as fol- 
lows: ‘A function of n variables Ly, Le, Ls, --:, Ln is said to be sym- 
metric in these variables if any interchange of the variables leaves the 
function the same. . . . Since any permutation of variables may be 
obtained by successive interchange of two variables, a necessary and 
sufficient condition that a function be a symmetric is that any inter- 
change of two variables leaves the function unaltered.’’ 4 

The nomenclature for a symmetric function has been established by 
prior usage.° 

A function f is called an SF if and only if f is equivalent to some 
function g where g is a symmetric function. Two functions are con- 
sidered equivalent when one may be obtained from the other by 
complementing some variables. 

When a function f is an SF, the variables of g (any symmetric 
function equivalent to f) and their complements are called COS’s. 
Such a pair defines an axis of symmetry uniquely specified by either 
member of the pair. Although any SF has at least two COS$’s, the func- 
tion is not symmetric in the same degree (the subscript in standard 
symmetric notation) about the two centers. In fact, if a function can 
be represented as m out of n about one COS, then it is (n — m) out of n 
about the complemented COS. 


* This circuit realization could have been as easily done with FET’s in MOS tech- 
nology. Each contact would be replaced with a single FET. 


™ The details of this are presented in the appendix. 
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Fig. 1—Simplified circuit for Born and Scidmore problem by Born and Scidmore 
technique. 


A function is an ASF if it 


(2) Contains one or more imbedded SF’s, i.e., the function can be 
expressed as the sum (oring) of two or more functions where 
one or more functions are SF’s. 

(zt) Is the complement of a function containing one or more im- 
bedded SF’s. 
(272) Can be expressed as the sum of two functions, one of which is 


the complement of a function containing one or more imbedded 
SF’s. 


Any function is an ASF if the limits are stretched far enough, but 
to benefit from symmetricity the number of imbedded SF’s should be 
small and the terms not included in the SF’s should be relatively few. 
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Fig. 2—Simplified circuit for Born and Scidmore problem by author’s technique. 
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Any function fits into one of three categories : 


(t) The function is not symmetrizable. 
(22) The function is symmetrizable about all points. 
(477) The function is symmetrizable about exactly two centers of 
symmetry.® 


The first category is not the subject of this paper. 

The second category consists of only four members, the ‘‘zero’’ 
function, the ‘‘one’’ function, and the odd and the even parity func- 
tions. The first two functions are trivial, while the parity functions 
are very specific functions that have been the subject of numerous 
papers.* These functions have exactly 2?-! minterms for a function of n 
variables, and no minterm differs from any other minterm in the state 
of an odd number of variables. 

The third category is of practical interest, being the class of SF’s 
that are not just trivially symmetric. It is also the large majority of 
SF’s with more than three variables (n > 38), since this category has 
27t1 — 4 members. 

A necessary and sufficient condition for a function to be symmetric 
about a specific COS may be expressed as follows: If one minterm 
matches the COS in exactly m out of the total of n variables, then ex- 
actly »C,' minterms must match that COS in exactly m out of n 
variables. This is a direct result of the definition of an SF. 

Thus, a simple test exists for ascertaining whcther a given function 
is an SF about a specific COS.* The problem then is determining what 
a COS must be if the function is an SF. Consequently, rather than 
exhaustively testing a function for symmetry, our technique first finds 
what the COS must be if the function is an SF and then verifies the 
actual symmetry about that point. 

A direct result of the theorem"! 


]= 
SE? (Liss, Lint, me “y L,)8 


is that a function f is an SF in Jy, Ly, L3, ---, L, only if a specific 
subset of the minterms is an SF in In, Le, Lz, +--+, Lm form S n. Thus, 


* Garner (Ref. 7) is one of many who have studied this function. 
mCn is the number of combinations of n things taken m at a time and equals 
n!/L(n — m)!n!). 
* In essence, this is the same test that McCluskey (Refs. 8 and 9), Marcus (Refs. 
5 and 10), and others have used. 
§ Si(Li, Le, Ls, +++, Ln) is standard symmetric notation for “symmetric a-out-of-n 
function of variables L1, Le, +++, Din.” 
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for any SF, a COS in n variables when all other variables are held con- 
stant is a subset of the total COS of the function. 

We have shown that a necessary condition for an SF is that the func- 
tion be symmetrizable in subsets of the variables. This, however, is 
not a sufficient condition even if the subsets encompass all variables, 
as we can show by an example. The eight-variable function 
f = Inleh3L483(Ls5, Le, L7, Lg) is symmetrizable in both the four most 
and the four least significant variables, but is not an SF in all eight 
variables. Even the fact that the function is symmetrizable in over- 
lapping subsets of variables is not a sufficient condition, as we can 
demonstrate by the seven-variable function f = S§(L1, Le, ---, L7) 
+ (Ly LeL3L',L;L,L7). Thus, the possibility exists that a function is 
not an SI’, even though a composite of the COS’s of subsets can be 
found. However, if COS’s are found for subsets of variables where the 
group of subsets includes all variables in the function, then a COS of 
the function (if it exists) must be the composite of the COS’s of the 
subsets. 

Since each COS has a mate (the point where all variables are the 
complements of the variables of the first COS), a COS of the func- 
tion must be a composite involving one or the other COS. The am- 
biguity as to which COS of each pair to select can be resolved by 
having each set of variables overlap another set in at least one vari- 
able. Thus, if there are n variables and each subset selected has & 
variables, then at least |(m — 1)/(& — 1)| subsets are required to 
find the COS.* 

S. H. Caldwell!:? has demonstrated a technique using Karnaugh maps 
for recognizing symmetry in functions of three or four variables. We 
review it here. 

Any single square on a Karnaugh map represents an SF of n out of 
n, Where n equals the number of variables. Thus, on a three-variable 
map, each square is an SF of three out of three [written S3(Z1, Le, Ls) | 
of some set of variables Ly, L2, and L;. As an example, the square a 
in Fig. 3 is A’BC’ and is the SF’ 83(A’BC’). 

On a Karnaugh map, only one variable changes state in going from 
one square to an adjacent square. Thus, the four adjacent squares 
match the center square in two out of three variables.’ Similarly, any 
squares that are two squares from the given square match it in one 
variable, etc. Thus, the squares labeled 2, 1, and 0 are S3, S?, and S§ of 








* 


| Notation for the least integer equal to or greater than the argument. 
* To see the pattern, it is necessary to extend the map-repeating columns and rows, 
as shown in Fig. 3. 
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Fig. 3—Symmetry in a three-variable Karnaugh map. 


square a, respectively. Note that the coordinates of square a are the 
variables of symmetry I, Le, Lz, etc. (i.e., the COS). 

Iixpanding this to four variables represents the situation shown in 
Fig. 4, where the 3’s represent $3 of square @, the 2’s represent S3 of 
square 8, etc. 

Since each SF of three or four variables yields a distinctive pattern, 
pattern recognition permits recognizing any SF of three or four vari- 





Fig. 4—Symmetry in a four-variable Karnaugh map. 


806 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1974 





4 ' = <4 ' 4 ’ 


Fig. 5—Sum of two fundamental symmetrics. 


ables, and the COS identifies the variables of symmetry. Note that 
Soi, D2, Ls, DL) -_ Si(Ly, Ln, Ls, I), ete. 

The identification of an SF which is the sum of more than one fun- 
damental SF in the same variables is depicted in Fig. 5. This method 
can be extended to more than four variables, as explained below. 

At this point, we depart from Caldwell’s technique. Note that 
Caldwell has presented a means of expanding his method to more 
than four variables (in theory, to any number). However, for 7 vari- 
ables, his technique requires the use of 2”-* Karnaugh maps and 2”-*—1 
applications of an expansion theorem." 

Since a technique exists for handling subsets of four or less variables, 
four can be substituted for & in the expression |(n — 1)/(k — 1)|, 
yielding | (mn — 1)/3]| as the minimum number of subsets required to 
find the COS of a function of n variables, if it exists. By appropriate 
selection of subsets, it is possible to require only the use of the mini- 
mum number of subsets. 


2.2 Technique 


Tirst, we list the terms in decimal notation* where unprimed and 
primed variables are represented by binary ones and zeros, respectively. 


“The use of decimal notation is not essential, but it speeds up the mapping and 
selection of terms which differ only in a specific set of successive variables. 
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Then, some four-variable Karnaugh maps are drawn. In general, one 
such map covering the lowest decimal numbers in the function forms a 
good starting point. Figure 6 presents a review of the decimal notation 
for Karnaugh maps with six variables. 

Once the Karnaugh map has been plotted, we can try to identify an 
SF on it. If the submap is quite full of ones, checking the remaining 
zeros for symmetry may be easier. The SF of S§124(A, B, C, D’) 
shown in Fig. 7 is easier to identify as the complement of 
S3(A, B, C, D’). 

If no COS is found on a map, the function is not an SF. If a COS 
is found, then at least one of the four variables involved plus not more 
than three new variables is plotted in a similar fashion.* The resulting 
terms are plotted on a new four-variable map, and the COS (if it 
exists) 1s found. 

This technique is repeated as many times as is necessary to account 
for all variables. Since the minimum number of maps is | (n — 1)/3], 
the minimum (and usual) number of maps for nine variables is three. 
Figure 8 tabulates the minimum number of extended four-variable 
Karnaugh maps required by this technique as compared with the 
Caldwell-Grea technique. Once the various COS’s are found, they 
are combined to form the COS of the whole function, complementing 
all the variables in any COS required to make the overlap variables 
match. This possibility exists since either of two COS’s could be 
selected—i.e., S3(ABCD) = S{(A'B’C’D’). 

Once the potential COS is found, each term of the function is com- 
pared with it to see if a complete SF is represented. 

The only restriction on the selection of a subset of minterms for 
which all variables are fixed except a specific four is that the subset 
must have at least one member and may not have either 8 or all 16 
members.' 

The reason for the eight-term restriction is that eight terms in four 
variables represent either the odd or even parity function, which are 


“A suggested technique for this is to divide the decimal value of each term by 2 
to include one new variable, by 4 to include two new variables, and by 8 to include 
three new variables. If insufficient terms are found by this technique, any number less 
than the divisor can be subtracted from each term’s value prior to the division. 
Division by 8 selects those terms ending in 000 and discards the last three terms. 
Subtraction of 2 followed by division by 8 selects the leftmost n — 3 bits of those 
terms ending in 010, ete. 

T Si) 23,4 (ABCDEFGHILJ) has one full subset for any four variables (i.e., the 
subset where all fixed variables are zero although all other subsets in those variables 
are not full). There are no terms greater than 640 in the above function, since all 
such terms require at least five ones. 
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Fig. 6—Decimal notation for Karnaugh maps. 
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THE O's REPRESENT S3 (ABCD’) 


THEREFORE THE 1's REPRESENT Sj , 5 ,(ABCD’) 


Fig. 7—Identifying a large symmetric. 
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NUMBER OF || NUMBER OF EXTENDED FOUR VARIABLE NUMBER OF APPLICATIONS 


VARIABLES KARNAUGH MAPS REQUIRED OF EXPANSION THEOREM 
4 1 1 0 0 
5 2 2 0 1 
6 2 4 0 3 
Z Z 8 0 7 
8 3 16 0 1 
9 ) 32 0 31 

10 3 64 0 63 
11 4 128 0 127 
12 4 256 0 255 
13 4 512 0 511 
14 5 1024 0 1023 
15 5 2048 0 2047 
16 5 4096 0 4095 
17 6 8192 0 8191 
18 6 16,384 0 16,383 
19 6 32,768 0 31,768 
20 t 65,536 0 65,535 


Fig. 8—Comparison of author’s technique and Caldwell-Grea technique. 


symmetrical about all possible minterms in four variables. In general, 
the selection of another set of values for the fixed variables results in 
a map without eight terms. * 


2.3 Illustrative example 


We test here a 10-variable function for symmetry. Trying all 2” 
possible combinations of terms or plotting the 64 four-variable Kar- 
naugh maps and mathematically combining them are unattractive. 
On the other hand, a 10-variable function is not unwieldy by this 
technique. 

Let us consider the function in variables A through J where 
F= ¥) (18, 21, 25, 28, 31, 37, 41, 44, 47, 49, 52, 55, 56, 59, 62, 93, 109, 
117, 121, 124, 127, 157, 173, 181, 185, 188, 191, 253, 285, 301, 309, 313, 


* S123 (ABCDEFGHTIJ) will exhibit apparent parity for only those subsets where 
all fixed variables have a value of zero. 
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316, 319, 381, 445, 541, 550, 557, 565, 569, 572, 575, 637, 701, 706, 829, 
834, 898, 960, 963, 966, 970, 978, 994). 

Figure 9 is a worksheet listing the terms of the function and the 
terms to be plotted on Karnaugh submaps. The appropriate Karnaugh 
submaps are shown in Fig. 10. 

For the subset in which A through F are zero, the only term is 13 
and thus it must be a COS of the four variables GHJJ (1101). Next, 
we find terms which fit the pattern 000X XX X000. Division by 8 yields 
one such term, i.e., 7, and so the COS in the four variables DEG 
must be O111. 

Next, we find a submap of the variables ABCD by selecting those 
terms ending in 111000 (i.e., subtracting 7 from the terms found by 
the first division and then dividing by 8). The only resulting term is 
0000. Thus, the COS of the whole function (if it exists) must be a 
composite of ABCD = 0000 and DEFG = 0111, and GHIJ = 1101, 
i.e., 0000111101. 

All that remains is to verify whether or not the function 1s an SF 
about this COS. Consequently, the COS is filled in at the head of the 
“MATCH” column on the worksheet and each individual term is 


A A 
n n/8 | (A-7)/8 TERM MATCH n n/8 | (A-7)/8 TERM MATCH 
13 301 
21 309 
25 313 
28 316 
31 319 
37 381 
41 445 
44 450 
47 541 
49 557 
52 565 
55 569 
56 7 0 572 
5g 575 
62 637 
93 701 
109 706 
117 829 
121 834 
124 898 
127 960 | 120 
157 963 
173 966 
181 970 
185 978 
188 994 
191 
253 
285 


Fig. 9—Worksheet for sample problem. 
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QOO0000XX XX SUBMAP 000X XX X000 SUBMAP 
17. 10 OO 01 


| | 
—+—4—4--+4--|--+-+4+- 
bot Fouwot tf 4 





COS = 1101 COS = 0111 


XXXX111000 SUBMAP 
11. 10 OO 0O1 





COS = 0000 


Fig. 10—Karnaugh submaps for sample problem. 


compared with this pattern and the number of matching variables 
recorded. Figure 11 presents the completed worksheet. Next, the 
number of matches is compared with the number required for an SF, 
1.€., mCn. For m = 8 and n = 10, »C,n = 10!/(8! X 2!) = 45 and for 
m= landn = 10,,C, = 10!/(1!9!) = 10. 

Since 45 terms match on eight variables and 10 terms match on one 
variable, the function is an SF, i.e., S!%(61).* | 

A technique has been presented that is simple to use manually and 
requires no extensive memorization of a routine such as required by 
the Marcus-McCluskey method.’-8 However, the real power of the 


“This is the decimal notation abbreviation for S!°s (A’B'C’D'EFGHI'J). 
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A MATCH A MATCH 


n n/8 | (A-7)/8 TERM 0000111101 n n/8 | (A-7)/8 TERM 0000111101 
13 0000001101 8 301 0100101101 8 
21 0000010101 8 309 0100110101 8 
25 0000011001 8 313 0100111001 8 
28 0000011100 8 316 0100111100 8 
31 0000011111 8 319 0100111111 8 
37 0000100101 8 381 0101111101 8 
41 0000101001 8 445 0110111101 8 
44 0000101100 8 450 0111000010 1 
47 0000101111 8 541 1000011101 8 
49 0000110001 8 557 1000101101 8 
52 0000110100 8 565 1000110101 8 
55 0000110111 8 569 1000111001 8 
56 7 0 0000111000 8 572 1000111100 8 
59 0000111011 8 575 1000111111 8 
62 0000111110 8 637 1001111101 8 
93 0001011101 8 701 1010111101 8 
109 0001101101 8 706 1011000010 1 
117 0001110101 8 829 1100111101 8 
121 0001111001 8 834 1101000010 1 
124 0001111100 8 898 11100000 10 1 
127 0001111111 8 960 | 120 1111000000 1 
157 0010011101 8 963 1111000011 1 
173 0010101101 8 966 1111000110 1 
181 0010110101 8 970 1111001010 1 
185 0010111001 8 978 1111010010 1 
188 0010111100 8 994 1111100010 1 
191 0010111111 8 
253 0011111101 8 
285 0100011101 8 


Fig. 11—Completed worksheet for sample problem. 


technique is that it can readily be extended to cases that are not pure 
SF’s, while the other methods cannot be so extended. The reason the 
new technique can be extended is because it depends on pattern recog- 
nition at which humans are adept. Thus, patterns can be discerned 
in spite of extraneous data, 1.e., “noise,’”’ while a technique which 
rigorously uses all data in a prescribed functional relationship cannot 
sort out the desired data. Obviously, an SF can be so obscured by 
‘noise’ as to be unrecognizable. However, the cases of the most value 
are those which have only limited obscuring terms. 

In the next section, some extensions of this technique are presented. 


Ill. EXTENSIONS OF THE TECHNIQUE 


Since the technique may be extended to cover a variety of situa- 
tions, we discuss some specific extensions here. The author has worked 
problems for each of the discussed extensions, and at least one of each 
has been presented in unpublished memoranda, while one example of 
a composite function is solved in the appendix. 
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For the purpose of this paper, a function is a multiple SF if it is the 
sum of two or more SF’s whose COS’s are not the same or com- 
plements of each other. Thus, the function f = S3(ABCDE’) 
+ §3(A’BCDE) is said to be a multiple SF. Since any minterm of a 
function of n variables is the SF Si (minterm), any function of m 
minterms can be represented as the sum of not more than m SF's. 
However, the useful cases are functions that are the sum of considerably 
less SF’s than minterms, such as the above example. 

An incompletely specified function is a function that contains at 
least one minterm for which transmission is neither required nor 
forbidden (don’t-care terms). 

A function is called an incomplete SF if all but a few minterms re- 
quired for an SF are present. 

An overly complete SF is a function composed of an SF plus addi- 
tional minterms. Any function can be forced to fit this mold, but the 
cases of interest are those in which almost all the minterms are included 
in the SF. 

Approximate SF’s are functions that are not truly SF’s but that 
differ only slightly from an SF. Incomplete SF’s are included in 
approximate SF’s, as are overly complete SE’s. However, approximate 
SF’s also cover functions that are the sum of incomplete SF’s and 
some additional minterms. 

Obviously, any function fits into the category of an approximate SF, 
but the useful cases are those whose deviation from an SF is slight. 
For example, the function f = S§8(ABC’DE'FG'’H) + AB’C’D' EF'GH' 
— ABC’D'EF'GH’ would fit this category. Functions like this can 
often be simply constructed by modifying the SF. 

A partial symmetric has been defined as a function which is an SF 
in some but not necessarily all variables.’ Thus, SF’s are a subset of 
partial symmetrics. However, usually a partial symmetric refers to one 
which is not an SF in all variables. 

Both multiple SF’s and overly complete SF’s are examples of func- 
tions containing imbedded SF(s). An incomplete SF 1s a function whose 
complement contains an imbedded SF, while approximate SF’s are 
the sum of two functions, one of which is the complement of a function 
containing an imbedded SF. 


3.1 Overly complete SF’s 


In many cases, a function can be expressed as an SF plus certain 
additional terms. That is to say, the function has an imbedded SF. The 
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eases of interest are those in which the number of additional terms is 
small, since a symmetric design can be used for most terms and then 
the additional terms can be treated individually. In some cases, part 
of the additional terms can be combined directly with the SF. 

The technique is basically the same as it 1s for a pure SF. However, 
more maps may be required because of the masking effect of potential 
extraneous minterms. An extra minterm on a map could result in an 
ambiguity as to where the COS is. Worse than this, a single extraneous 
term could be the only term which appears on a specific submap. This 
could lead to the selection of a false COS for the entire function and an 
indication that the function did not contain an SF. 

Because of this, it is best to select two submaps in the same four 
variables and get a concurrence as to the COS of the four variables. It 
is, of course, possible to find in some functions pairs of submaps In the 
same variables which do not concur on a single COS. Then a third 
submap is required in those same four variables. In the author’s ex- 
perience, this seldom occurs. 

Theoretically, it is possible to have to continue making more sub- 
maps without actually determining what the COS of the four vari- 
ables must be. However, the cases of the most value are those cases 
that require the least maps. In practice, if concurrence on a COS can- 
not be found by the use of four submaps in the same four variables, no 
very useful SF exists within the function. Also, if any one submap 
has more than three terms, it is usually possible to determine the COS 
of the four variables. 

When dealing with an SF, we stated that a submap with eight terms 
represents one of the two parity functions. The parity function is 
recognized by the fact that, on a Karnaugh map, no two adjacent 
squares have the same value (zero or one). If eight terms appear on 
the submap and this relationship is not true, then the COS can usually 
be readily determined. 

As before, once the various COS8’s are found, they are combined to 
form the point which must be the COS if the function is primarily 
composed of an SF. Then each term in the original function is com- 
pared with it to determine if all the terms of an SF are represented and 
which, if any, additional terms are included in the function. 


3.2 Multiple SF’s 


It is usually possible, if there are two or more SF’s imbedded in a 
function, to find each of them. The technique is based on the principle 
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used in finding an SF. Again, submaps are selected for the sets of four 
variables, using more than one submap in the same four variables if 
necessary. In practice, there is often less confusion when the additional 
terms do form an SI’. Usually it is possible to see more than one COS 
in a single submap unless the two SI'’s both have the same COS in the 
chosen variables. However, since it is possible to have two SF’s which 
have the same COS in four variables although not in all the function 
variables, the appearance of only one COS in a submap does not rule 
out the possibility of having more than one SI’. Also, since in some sub- 
maps no terms of one SI* may appear, if a function is found to have 
some isolated terms, these terms should be investigated to see if they 
(plus perhaps some terms in the discovered SI) form another SI’. 

At any point in the analysis of the function that more than one SF 
appears to be found on the same submap, the remainder of the analysis 
should be done on the basis of a supposed multiple SF. When two (or 
more) COS8’s are found that are not identical or mates, each COS is 
used as a basis for determining the minterms to be used in a submap in 
the next selected set of variables. In general, selection of two new 
variables and two old variables works better when a multiple SF is 
suspected, since the new set of variables must contain enough informa- 
tion to make possible the selection of the appropriate COS. As long 
as the variables to be held fixed while new submaps are defined differ 
in at least one position, different submaps can be defined for finding the 
next COS, i.e., the submaps 0OX XX X11 and 0OOX XX X01 are different 
and usually will give COS’s that can be readily correlated with the 
COS8’s found in the lowest-order submaps (perhaps 1011 and 0001). 


3.3 Incomplete SF’s 


An incomplete SF is a function that lacks a few specific terms of 
being an SF. Thus, it is the complement of a function with an 1m- 
bedded SF. Such a function can be expressed as one of the three fol- 
lowing functions. 


(z) A modified symmetric circuit. 

(21) A symmetric circuit ANDed with another circuit which blocks 
those terms supplied in the symmetric circuit that are not part 
of the total function. 

(zz) A combination of the above. 


In the first case, such a modified symmetric circuit would be very 
efficient. In the second and third cases, the value of finding the SF 
decreases as the complexity of the blocking circuit increases. 


816 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1974 


The technique for finding an incomplete SF is similar to that used in 
previous examples. The major point of difference is that some terms 
needed to complete the SF are not in the original function, and hence 
they must be added to the function to form an SF and then blocked 
by circuit changes. 

Finding the COS may require addition of certain terms to the 
function. Once the proposed COS is found, additional terms may be 
required to complete the SF. 

All terms added to find the COS or to complete the SF must be 
recorded so the effect of these terms can be deleted in the realization 
of the function. 


3.4 Incompletely specified functions 


An incompletely specified function is one that has at least one 
minterm for which the function is undefined. The use of some terms 
for which the function is undefined in conjunction with those terms for 
which the function is defined often permits simpler circuit configura- 
tions than could be achieved if all these extra terms are ignored (tacitly 
forced to a definition of no transmission states). The literature abounds 
with references to use of these terms in conjunction with Karnaugh 
maps. Although the use of these terms to aid in completing an SF has 
been previously ignored by most authors, some work has been done by 
Arnold and Lawler.® 

Since the method described here depends on pattern recognition, 
the use of some terms to complete map patterns is straightforward. 
Once some (or none) of the undefined minterms have been used to 
ascertain what the COS of the function must be, then all terms for 
which transmission is required plus all terms for which the function 
is undefined are compared with the COS. Those undefined terms which 
have the same order of symmetry as the required terms are then used 
to test for the symmetricity of the function and, if sufficient terms are 
found, these terms are used to transform the function into an SF. 


3.5 Approximate SF’s 


An approximate SF is a function that differs only slightly from a 
true SF. The cases of greatest interest are those that have a few terms 
not in the SF and are missing some terms needed to complete the SF. 
Thus, we can think of this as a combination of an overly complete SF 
and an incomplete SF. The most interesting problem solutions in these 
functions occur where the extra and the missing terms can be paired 
to result in only a minor modification to the true SF. 
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Here, as before, the basic attack is to find a COS about which the 
majority of the terms must be symmetric if the function is to approxi- 
mate an SF. Since both missing and extra terms may occur, it may be 
necessary to ignore some terms in attempting to find a COS and it 
may also be necessary to temporarily add to the function certain other 
terms to complete the symmetry in some subset of variables. Once a 
proposed COS is found, certain terms may not be symmetric about 
that COS, and there may not be enough terms symmetric about the 
proposed COS to complete the SF. 


3.6 Composite function 


A composite function is a function that may be an incompletely 
specified function and that may contain one or more imbedded, in- 
complete, or complete SF’s. In other words, the function can be almost 
any function. 

The procedure is basically the same as has already been used. How- 
ever, because of the possible complexity of the function, some feature 
may be obscured and, even though a COS for part of the terms is found, 
it may be desirable to repeat the process, treating all terms in the first- 
found incomplete or complete SF as don’t-care terms for subsequent 
analysis. This procedure can be repeated until the work entailed is not 
warranted by the number of remaining terms not identified with an 
SF. The Born-Scidmore® problem considered in the appendix is an ex~ 
ample of a composite function. 


IV. CONCLUSIONS 


A technique has been presented here for isolating SF’s (complete or 
incomplete) in the presence of other SF’s or other terms, or in in- 
completely specified functions. 

Since the method depends basically on pattern recognition, at which 
humans excel, as opposed to value manipulation, at which they do 
not excel, this technique is especially useful for manual use. However, 
since the number of patterns is quite restricted, such a technique could 
be implemented by a computer program. 

Since symmetrics are powerful tools in implementing functions, the 
recognition of SF’s within functions can greatly reduce the work in the 
synthesis of functions. 
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APPENDIX 


In Reference 3, Born and Scidmore have transformed a partial 
symmetric in six variables into an SF (pure symmetric) in nine 
variables. Specifically, the function they examined was F = }~ (1, 2, 3, 
4,5, 6, 9, 10, 11, 12, 18, 14, 17, 18, 19, 20, 21, 22, 24, 25, 26, 28, 33, 34, 
35, 36, 37, 38, 40, 41, 42, 44, 48, 49, 50, 52, 56, 57, 58, 60). The minimum 
SF they presented was S33,5 (A, B, C, D, D, HE, E, F, F). This SF 
can be represented as shown in Fig. 12 and, as it is shown, requires 
34 transfers plus 2 break contacts with one relay requiring 10 transfers 
and 2 breaks. 

However, this circuit can be simplified by standard techniques. The 
resulting simplified circuit is shown in Fig. 13, which uses 28 transfers 
and 2 breaks with the largest single-contact load being 5 transfers 
and 1 break. 

This problem can also be solved by the author’s technique. Figure 
14 is the worksheet and Fig. 15 the first set of Karnaugh submaps. 
The 0OXXXX submap shows two axes of symmetry, one with one 
center at 0111 and the other with one center at 1111. Subtracting 3 and 
dividing by 4 yields four terms for the XXX X11 submap. These four 
terms almost identify a center of symmetry at 0000 (or 1111). Either 
one term is missing (0001) which would have come from a minterm 7 
or the term 0000 is actually that term shifted out of place by a modifi- 
cation of the symmetric. Either case argues for the use of this 1111 as 
a center of symmetry. When this is combined with the centers of sym- 
metry in OOX XXX, one match is found, namely, 111111. Thus, the 
terms are compared to this center of symmetry and 15 terms are found 
to match it in two variables. Since this is the value of Cs, the func- 
tion contains S$(63). Next, the remaining terms are plotted on the 


RELAY CONTACT LOAD 





mmOoOQdOBwY 


TOTAL 


AD BS Ae D oD «E E FF 
Fig. 12—Symmetric for 83,345 (A, B, C, D, D, E, E, F, F). 
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RELAY CONTACT LOAD 


A 4T 
5T 
5T,1B 
57T,1B 
5T 
4T 


TOTAL 28T,2B 


AMONG 





A E F AB Cc 80D E F 
Fig. 13—Simplified symmetric for S3,3,45 (A, B, C, D, D, E, E, F, F). 


MATCH | MATCH TERMS MATCHING 

N |(N-3)/4] TERM | 111111 1110111 |NEITHER IN TWO VARIABLES 
1 000001 1 2 

2 000010 1 2 

3 0 000011 2 3 

4 000100 1 2 

5 000101 2 3 

6 000110 2 3 

9 001001 2 1 

10 001010 2 1 

11 2 001011 3 2 

12 001100 2 1 

13 001101 3 2 

14 001110 3 2 

17 010001 2 3 

18 010010 2 3 

19 4 010011 3 4 y 
20 010100 2 3 
21 010101 3 4 J 
22 010110 3 4 J 
24 011000 2 1 
25 011001 3 2 
26 011010 3 2 
28 011100 3 2 
33 100001 2 3 
34 100010 2 3 
35 8 100011 3 4 y 
36 100100 2 3 
37 100101 3 4 / 
38 100110 3 4 / 
40 101000 2 1 
41 101001 3 2 
42 101010 3 2 
44 101100 3 2 
48 110000 2 3 
49 110001 3 4 y 
50 110010 3 4 J 
52 110100 3 4 / 
56 111000 3 2 
57 111001 4 3 J 
58 111010 4 3 J 
60 111100 4 3 / 

MISSING TERMS 

16 010000 1 2 
32 100000 1 2 


Fig. 14—Worksheet for Born and Scidmore problem. 
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OOXXXX SUBMAP XXXX11 SUBMAP 
00 6~O1 6©611:6©61006000=CO01 SS 11—Ss«*210 





COS = 0111,1111 COS = 0000,1111 
Fig. 15—Karnaugh submaps of Born and Scidmore problem. 


submap OOXXXX (Fig. 16), and a single axis of symmetry (center 
0111) is found.* The same last three digits are here as before, so the 
same terms appear on the X XX X11 submap except for the terms which 
came from S8(63). 

If the three terms on the XXXX11 submap define a center of 
symmetry with a missing term, then this center of symmetry must be 
the same as before. Part of the minterms of S$(63), S?(63), and S?(63) 
appear, but none in sufficient quantity. Therefore, we assume that 
two of the terms on the XXX X11 submap are from extra terms and 
try each of the three points as a center of symmetry. Only one of these 
(1101) can be combined with 0111 so that a proposed center of sym- 
metry is 110111. The terms are compared with it and 13 terms are 
found to match 110111 on two variables. Since 15 terms are required, 
the comparison is made for the missing terms. With experience, this 
can be done very readily. However, for completeness, the terms of 
S3(55) are tabulated in Fig. 17. The missing terms are 16 and 32. Thus, 
the function contains the function $§(55)-16-32. 

Looking at the 12 remaining terms we note that three of them (57, 
58, 60) can be accounted for as ABCS?(DEF). The remaining nine 
terms all have the third variable in the zero state and so can be ex- 
pressed as C’ f(ABDEFP). 

Repeating the technique on the nine-term, five-variable function 
results in the OXXXX, IXXXX, XXXX0O, and XXXXI1 submaps 
shown in Fig. 18. The first two submaps concur on a 1111 COS and the 


“ Note that the terms already identified appear as ‘‘don’t care’s’’ (D) in the sub- 
maps even though, in this case, they do not result in simplifying the remaining im- 
bedded functions. 
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OOXXXX SUBMAP XXXXI1 SUBMAP 
0o0 O01 #11 #=10 OO O11 1%t 10 00 O01 #141 #=+10 OO OT 11 10 





COS = 0111 COS = 0010/1101,0100/1011,1000/0111 


Fig. 16—Karnaugh submaps of remainder of function. 


latter two submaps also result in a 1111 COS. As a result, a composite 
COS of 11111 is tried. Nine of the 10 needed terms for S$3(31) are 
present with only the term 00111 missing. Thus, the original function 
has been broken up into the sum of four smaller functions, namely: 


(i) S§(ABCDEP). 

(ii) S§(ABC'DEF) — A'BC'D'E'F' — AB‘C'D'E'F’. 
(iii) ABCS3(DEF). 

(iv) C'S8(ABDEF) — A'B'C’DEF. 


INCLUDED IN FUNCTION 


111000 
100000 
101100 
101010 
101001 
010000 
011100 
011010 
011001 
000100 
000010 
000001 
001110 
001101 
001011 


v 
v 
v 
v 
v 
v 
J 
v 
v 





Fig. 17—-Terms of symmetric S§ (55). 
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IXXXX SUBMAP 


OXXXX SUBMAP 








COS = 1111 


COS = 1711 


XXXX1 SUBMAP 


XXX XO SUBMAP 


10 


11 


00 01 


10 


01 11 #10 


10 00 


1 | 


Q 


1 4 








COS = 1111 


COS = 1111 


= 11111 FOR ABDEF 


COS 
MISSING TERM = 00111 


Fig. 18—Karnaugh submaps for final nine terms of function. 


Fig. 19—Circuit for S§ (63). 
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- BLOCK E'D'F'BA‘C’ 
i INSERT B BREAK HERE 





B 
F E D B A.C 
/ 
BLOCK E’D’F’B’AC’ ~ 
Fig. 20—Construction of circuit for S§(55)-16-32. 
F E D B A Cc 
Fig. 21—Symmetric for ABCS? (DEF). 
F E D B A Cc 


Fig. 22—Symmetric for C’S}(ABDEF)—A'B’C' DEF. 
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RELAY LOAD 








AMOoUOWDW>DYD 
Ce | 





TOTAL 117,18 | 





F a D B A Cc 


Fig. 23-—Final circuit for Born and Scidmore problem. 


These represent 


(7) An SF. 
(zt) An incomplete SF. 
(227) A minterm in three variables ANpbed with an SF in the re- 
maining three variables. 
(zv) Asingle variable ANped with an incomplete SF in the remaining 
five variables. 


Figures 19 through 22 portray a relay configuration for each of these 
smaller functions, while Fig. 23 is the resulting circuit configuration 
when the preceding four circuits are combined. This resulting circuit 
uses 11 transfers and 1 break-contact, or less than half of the number 
required for the circuit resulting from the Born and Scidmore technique. 
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Finite-Element Method for the Determination 
of Thermal Stresses in Anisotropic Solids 
of Revolution 


By S. KAUFMAN 
(Manuscript received October 19, 1973) 


This paper discusses stresses and deflections in anisotropic solid struc- 
tures of revolution. It presents two methods based on finite-element tech- 
niques: one, a solid-of-revolution method in which material properties, 
applied forces, and temperatures are independent of angle, and two, a 
long-cylinder method in which these properties are independent of the longi- 
tudinal coordinate. These methods are postulated on uniform stress fields 
within the element, rather than on the usual functional displacement 
description within the element. A Fortran program has been written for 
both these methods, and ample test problems are presented to validate the 
methods. An application 1s presented for thermal stresses induced during 
the post-growth cooling stage of Czochralski-grown lithium tantalate 
crystals. 


1. INTRODUCTION 


This paper considers two finite-element networks. The first network 
consists of triangular annuli forming a solid of revolution of any 
arbitrary cross section in the radial-longitudinal plane. The network 
has 27 symmetry with regard to material properties, external loadings, 
and temperatures. The second network consists of trapezoidal and 
triangular elements in the plane of the circle forming a right circular 
(actually, a polygon cross section) cylinder long in the longitudinal 
direction. The restriction of 27 symmetry is lifted for the second 
network and is replaced with the restriction that material properties, 
external loadings, and temperatures are independent of the longitudinal 
direction of the cylinder. 

Both methods were programmed on the IBM 370 computer. For 
the first method, two test plane strain problems are presented: one, a 
hollow cylinder subjected to a negative radial pressure and two, a 
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solid cylinder subjected to a linear radial temperature gradient. Com- 
parisons of stresses and deflections with known plane strain solutions 
are excellent in both cases. In the second method, the solution of a 
long cylinder subjected to a linear radial thermal gradient is presented. 
Comparison with a known theoretical solution again was excellent. 

Results are presented for thermal stresses induced in a lithium tan- 
talate crystal during the post-growth cooling stage. Lithium tantalate 
is in crystal class C3,. By aligning the trigonal axis of the crystal with 
the longitudinal axis of the cylinder, the problem can be analyzed by 
both methods. Comparison of thermal stresses obtained from the two 
methods was very good. 


Hl. SOLID-OF-REVOLUTION METHOD 


The basic element for this method is the triangular annulus shown 
in Fig. 1. The element, defined in the radial-longitudinal plane (r — 2), 
has 27 symmetry. A network of these elements will comprise any 
desired solid of revolution, whether it is solid or hollow, cylindrical or 
conical, or any combination thereof. Material properties, temperature 
distributions, and external forces must be independent of angle. The 
material properties are then limited to an orthotropic system with 
isotropic properties in the plane of the circle. This limitation is lifted 
in the case of the long cylinder method presented in the next section. 


RADIAL 
v7 AXIS 
“ 





A="L i x Lig A=|a| 


Ai = V/L, {HALF ALTITUDE AREA (k TO ij)] 


ap=ep-e, b= ey: ey 


Fig. 1—Properties of triangular annulus. 
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Solid triangular elements have been used successfully by other 
authors;!:? however, a different approach is presented here. Rather 
than to assume a displacement function,!:? a simpler and more direct 
approach is to postulate uniform stress fields in the annulus. This 
method was applied successfully in the mid-60’s by David B. Hall for 
the triangular membrane, but unfortunately his work has not been 
published. Hall’s triangular membrane is extended here to a three- 
dimensional model by incorporating a uniform hoop stress in the 
annulus. 

The initial strategy is to relate equilibrium between forces at the 
three points of the triangle 2, 7, k, and stress fields parallel to the sides 
of the triangle o;;, oj, ox:, and the hoop stress o». Before we do this, 
let us first compute a few properties of the triangular annulus. These 
properties include the area of the triangle A, length of each side L,;, 
etc., volume V, areas A,;, etc. associated with the stress fields o;;, 
etc., and direction cosines a;;, b;;, etc. These properties are presented 
in Fig. 1. 

In matrix notation, equilibrium for annulus n between grid point 
forces and the four stress fields described above is given as follows. 


{Fa} = [LB] {oo}, (1) 
where 
[Fal = (PiP iP iP PaeP ec} 
oo} = {oijojonsoo}, 
and 
— Aisi; 0 Andis 20A/3 
—A ij; 0 A riQes 0 
_ A ijba; ae A jxb 5% 0 27rA/3 
Le 7 AijQij — A jpO;x 0 0. 
0 A jxb jx — Axi; 27A/3 
0 A jrQjx —AnriQks 0 


The r and z subscripts attached to the forces F;,, Fiz, ete. refer to 
the radial and longitudinal directions, respectively. 

The skew stress field o:;, oj, o%¢ 18 merely an artifice to allow us to 
readily establish the equilibrium relationship. We are actually inter- 
ested in the orthogonal field o., o,, 7,2, aS Shown in Fig. 2. The relation- 
ship between the two fields including the hoop stress is given below. 


{on} = LD ]{oo}, (2) 
where 


{ Tn} = { 720 rT 206} 
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res oe 
€a 
eK = 
e ij 
Fig. 2—Stress fields. 
and 
Oj; je Cz; 0 


bj; Di bi 0 
A:jb5; 540 5% AniDes 0 
0 0 0 l 


= 


After inverting eq. (2) and substituting it into eq. (1), the following 
relationship is obtained. 


{F,} = Pelee: (3) 
where 


[F] = (B]LDS". 


A conjugate relationship to eq. (3) by an application of the principle 
of virtual work can be stated as 


fe} = GLP), (4) 


where the strains are 
{én} = {€c€rVre€o}, 
the deflections conjugate to {F,} are 
(Ot "Ug 40 0 ee U eet: 


and the superscript ¢ denotes a transposition. 
The stress-strain relationship, including the thermal strains fad7, 


IS given as 
{on} = [C1( {en} - | fear), (5) 
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where [C] is a symmetric 4 X 4 stress-strain matrix such that 


Ch = C12, Cy = C'o4, Cig = Co3 = Cg3 = 0, 


| fear ~ | feud? far 0 faa}, a, = a0. 


Substituting eq. (5) into eq. (4) obtains the stresses in terms of the 
unknown deflections and known thermal strains as 


lon} = CO GLFIU.}—| faar| ) 6) 


Substituting eq. (6) into eq. (3) obtains contributions to the network 
stiffness matrix and thermal load vector for annulus 7 or 


{Fn} = [Kal{Un} — {En}, (7) 


and 


where 


[Kp] = SLFILCILFY 


is a symmetric 6 X 6 stiffness matrix and 


(B,) = CFIC] f aa? 


is a6 X 1 thermal load vector. 
Annuli contributions to the network stiffness matrix and thermal 
load vector are additive, and hence eq. (7) can be written for the net- 


work as 
{F} = [LK ]{U} — {#}, (7a) 
where 


LA] = 2 [K,] and {Bh} = 2, (Ba, 


and where the number of equations equals twice the number of points 
of the network (one longitudinal and one radial degree of freedom per 
point). 

The [K ] matrix in eq. (7a) is singular, as there is no constraint to 
prevent rigid body motion in the longitudinal direction. Rigid-body 
motion in the radial direction is prevented by the hoop stress field, as 
can be seen in eq. (1). (Note that in Eq. (1) os» is the only non-self- 
equilibrating stress field.) In addition to at least one reference longi- 
tudinal constraint, radial constraints must be provided for solid 
cylinders at points of zero radius to prevent radial (or hoop) motions 
at these singularities. The degrees of freedom can now be partitioned 
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into an unconstrained set denoted by ‘“‘a’’ and a constrained set de- 
noted by ‘b.” This partitioning can be represented schematically as 


PF, _ Forse U = De : (8) 
Fy| KK op UO, = 0 Ey 
From the upper equation of (8), the previously unknown deflections 


can now be obtained in terms of known external forces and thermal 
loads, or 


{Ua} = [Kaa l"({ Fa} + {£a}). (9) 


From the lower equation of (8) and with the help of eq. (9), the 
forces of constraint, if desired, can be obtained as 


{Fo} = (Kas) Kaal {Pa} — {Fo} — ({2o} 
— [Kas ][Kaal*{Ha}), (9a) 


where {F,} represents external forces applied at the constraint points. 


Y, en 





x, €, 
ne : 6 eet es 
a = tan Mtge th) sin AP Hrg—1,)1 df BB b= @-e 
@43 =[-cos(a—p)e,, sin (a— 9) ] ag = €@- e, bg = eg -&, 
€o4 = [-cos(a+O)e,, -sin(a+é)e, ] 4437 €,4°4, biz = €)3: & 
_»A80 Yi : fo en 
22 =[sin? = (rg+ 14)? + (rgtr)2 17 ang = On4° @, bog = 54° ep 
See eee : Ag AO 
e, = (cos #e., sine, ) A= (tg—t;) (rot r)) sin =— cos =5— (AREA) 
ey = (—sin@ e,, cos 6 e) 2c = (rg—F,) cos AO f= A/4c 


Fig. 3—Properties of the trapezoidal element. 
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A=K, 2 sine cosAP (AREA) 


2g=r, sinAé 


A@é Aé ] 


ey = [cos =e e,, sin (9 — “5 ) e ep 


ao = €10°@, bya = C12 & 
€3 = (—sin@ e,, cos 6 e,) 
4937 €53°@, bgg= C93: &b 


Aé 


cor = [- cos(@ +=) e,, —sin(@ AO 6,) 
31 @3,°€, 3, = €31,° & 


Fig. 4—Properties of the triangular element. 


lll. LONG-CYLINDER METHOD 


In this section, the 27 symmetry requirements as to material proper- 
ties, external loadings, and temperature are lifted and a two-dimen- 
sional model is constructed in the plane of the circle. The cylinder is 
considered long in the z-(longitudinal) direction and material proper- 
ties, external loadings, and temperatures are assumed independent of z. 
Shear stresses r,z and 72; are assumed to be zero. The most general 
stress-strain relationship considered is the 4 X 4 submatrix bounded 
by the dashed lines of eq. (18) (appendix) with the additional proviso 
of no coupling between the stresses oz, oy, oz, Tzy and the strains yyz, 
Vex. Tentatively, we let «, vanish during the initial phase of the analysis. 
This restriction is subsequently removed in a manner similar to that 
described by Timoshenko? for long isotropic cylinders. __ 

The basic building blocks of the long cylinder is the: trapezoidal 
element (Fig. 3) and the isosceles triangular element (Fig. 4, solid 
cylinders only). The trapezoidal element is subjected to three constant 
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stress fields: a radial stress o,, a tangential stress a», and a shear stress 
vr. The triangular element is likewise subjected to three uniform 
stress fields, each parallel to a side of the triangle. Equilibrium between 
forces collected at the apices of the elements and three uniform stress 
fields are obtained as shown below (see Figs. 3 and 4 for the properties 
of the elements). 


Trapezoid: 
| OP fa, —Cds lay3 
Py fo, — che lby3 
Fi. fa, Cae =< laos - 
F'y2 _ fo, cb =< lbos J 
Pus f= fa, cag = las ee (10) 
Fy = fob, cb —* lbi3 Hat 
Poa — fa, — Cae larg 
Fy rs fo, —cbe ldo4 
Triangle : 
Pay —  JaQyj2 0 ga31 
Py —gbre 0 gbs1 o 
Py. = gaye —1q/2 0 
Fy» 7 goie —1q/2 0 i (11) 
F.3 0 a/2 — 9Q31 ast 
F'y3 0 ta/2 —gbs1 


The relationship between the orthogonal stress field (¢,c97,-«) and 
the skew stress field (012 o23 o31) is given below for the triangle. 


Or 
{os} = 406 
Tré 
cos? A@/2 0 cos? A@/2 712 
= sin? A@/2 ] sin? A@/2 O23 7° (11a) 
—gsin A@/2 cos A6/2 O sin A@/2 cos A@/2| | oa1 


After inverting eq. (11a) and substituting it into eq. (11), equilibrium 
is established between the orthogonal stress field and the forces (de- 
noted as {F’,}) at the apices of the element or the points of the network. 
Denote this relationship as: 


{Fn} = LHal{oo}. (12) 
Similarly, for the trapezoid, eq. (12) will be the shorthand matrix 
notation for eq. (10). 
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Conjugate to the forces {F,} are deflections {U,} and therefore the 
conjugate relationship to eq. (12) can be constructed as 


{es} = 1/ALH,]{{ Un}, (13) 


where A is the volume (area times unit thickness) of the trapezoid or 
triangle. 

Upon substituting eq. (13) into eq. (26), the stresses in the element 
in terms of the unknown deflections at the apices and the known thermal 
strains are obtained. 


fos} = [0W(1/ALH (On) - CBs]| i aa}. (14) 


Upon substituting eq. (14) into eq. (12), contributions to the net- 
work stiffness matrix and thermal load vector are obtained for the 
element 

LE at = Leet C al _ {En}, (15) 


where 


[A,] = 1/A(H, ILC ILA, J 


is a symmetric 8 X 8 or 6 X 6 stiffness matrix and 


[B.) = (HaICCILES]| [oat 


is an 8 X 1 or 6 X 1 thermal load vector. 

Contributions of each element to the network stiffness matrix and 
thermal load vector are additive, and hence eq. (7a) can represent the 
force balance at all the points of the long-cylinder network as well as 
those of the solid-of-revolution network. 

Noting the material restrictions outlined in the beginning of this 
section, a network of triangles and trapezoids need only occupy one- 
quarter of a circle with planes of symmetry at @ = 0 and 7/2. Hence, 
tangential deflections must vanish at all points along these planes. The 
previously unknown deflections and constraint forces can now be 
solved in terms of known external forces and thermal loads in the same 
manner as was accomplished in the preceding section [see eqs. (8), 
(9), and (9a) ]. 

We are not quite finished. Recall that we have let «, = 0 throughout 
the cylinder, resulting in an axial stress a, [eq. (23) ] applied to the 
ends of the cylinder. If we superimpose a uniform axial stress (c,) 
such that the resulting force on the ends of the cylinder is zero, the 
self-equilibrating distribution remaining on the ends will, by St. 
Venant’s principle,* give rise only to local effects at the ends. The uni- 
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form axial stress correction to be added to eq. (23) is given as 


Cz Leq. aes. 
(oz) 7 (16) 


faa 


This correction results in added radial and tangential strains ob- 
tained from eq. (22). 


(e,(0)) = E13(0){oz) = (Ez. cos?@ + Hy, sin?8)(oz) (16a) 
(e9(0)) = Eo3(0){o2) = (zz sin26 + Ey, cos?) (cz). (16b) 
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Fig. 5—Thermal stresses plane strain problem. 
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The resulting correction to the radial and tangential deflections are 
next obtained. 


(U,(r, 9)) = r(e,(0)) = r(Ezy cos?6 + E,, sin?6) (oz) (17) 


(Tole, )) = f° (rveal0)) — (Uelr, 6)))a0 


; i " (y5(8) — E4s(8))(o=)d6 
= rsiné cosé (Hy, — Hz). (17a) 


Note that the correction (U,(r, 6)) vanishes at 6 = 0 and @ = 7/2, 
agreeing with the stipulated boundary conditions stated previously. 


IV. COMPARISONS WITH THEORETICAL SOLUTIONS 


Two plane strain problems are presented for the solid-of-revolution 
method: one, a linear radial thermal gradient applied to an isotropic 
solid cylinder (Fig. 5) and two, a negative pressure applied to the out- 
side circumference of a hollow cylinder (Fig. 6). Comparisons were 
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Fig. 6—Plane strain solution—thick-walled cylinder. 
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made of the finite element computer results with known theoretical? 
solutions. The agreement was excellent. 

A thermal problem is presented for the long-cylinder method in 
which an isotropic cylinder is subjected to a constant temperature 
plus a radial thermal gradient (Fig. 7). Comparisons were made of 
the finite-element computer results with a known theoretical solution. 
The known solution on page 410 of Ref. 3 was found to be in error. In 
the expression for radial displacement (U), insert (1 — 3v)/(1 — v) 
for (1 — 2v) and in the expression for c,, replace 2r with 2. After the 
above corrections were made to the theoretical solution, excellent com- 
parisons resulted as shown in Fig. 7. 


V. THERMAL STRESSES IN LITHIUM TANTALATE CRYSTALS 


Thermal stresses were computed for a lithium tantalate crystal, 
class C’3,, during the post-growth cooling stage. The crystal was ana- 
lyzed by the two methods presented in this paper on the IBM 370 
computer. For the solid-of-revolution method (Fig. 8), the model con- 
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Fig. 7—Thermal stresses—long cylinder. 
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Fig. 8—Network of triangular annuli—lithium tantalate crystal—solid-of- 
revolution model. 


sists of a right cylinder of radius a = 0.4875 inch and length 6 = 3 
inches plus a cone of length 0.875 inch attached to the far (cold) end. 
Radial constraints are provided at all zero radius points. For the long- 
cylinder method (Fig. 9), the model comprises one-quarter of a circle 
with tangential constraints at all points along the two planes of sym- 
metry (6 = 0 and 7/2). 





—— — 0.4375 — INCH RADIUS — — 
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Fig. 9—Network of triangles and trapezoids—lithium tantalate crystal—long- 
cylinder model. 
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Elastic properties’ of lithium tantalate are shown in Fig. 10. Axis 
1 is perpendicular to the vertical mirror plane and axis 3 corresponds 
to the trigonal axis of the crystal. If we strike out the small shear 
coupling terms, c14, in Fig. 10, then we have met the stress-strain 
relationship [eqs. (5) and (18) ] requirements of this paper. The 
coefficients of thermal expansion will be assumed constant in the 1, 2 
directions and piecewise linear along the trigonal axis 3. These co- 
efficients® are given below. 


ay = a. = 21.9 XK 10-8/°C 
as = 8.35 X 10-8/°C 625°C and higher 
— 3.1 K 107§/°C 500°C — 625°C 
0 400°C — 500°C 
1.16 X 10-§/°C = room temperature — 400°C. 


The crystal is assumed to be in a stress-free condition at an elevated 
temperature distribution 7'1, and elastic material behavior is assumed 
linear from the initial state of strain at temperature distribution 7'1 
to the final state of stress and strain at room temperature (26°C). 
This assumption rules out plasticity and stress relaxation. It is felt 
that ignoring these nonlinear effects plus ignoring the initial state of 
stress at 7'1 does not distort the qualitative picture of the stress 
pattern after the cooling-down process. The elevated temperature 
distribution in degrees centigrade is given below. 


T1 = 1500 — 100 (r/a)? — 500 z/b, for solid of revolution 
1500 — 100 (r/a)?, for the long cylinder. 


I 


The addition of the longitudinal gradient for the solid-of-revolution 
method was found to induce negligible stress in the crystal. 

For the solid-of-revolution method (z-growth), material axes 1, 2, 
and 3 correspond to r, 6, and z, respectively. For the long-cylinder 
method, two cases were investigated: one, z-growth where material 
axes 1, 2, and 38 correspond to 2, y, and z, respectively, and two, y- 
growth where material axes 1, 2, and 3 correspond to z, xz, and y, 


C33 = 2.75 x 10"! N/m? 


oy Gia) Cio at 0 0 E, Cy, 72.33 
O95 Cro O44 84971 0 60 E5 C14 =0.11 
73 — | 43 13 33 9 0 OD Ea Cy5 = 0.47 
es Dee Oncga: 0. 20 ¥93 C,3 =0.80 
731 0 0 OQ 0 C44 Or Y31 Ca, = 0.94 
Ty 0 0 0 0 SE CE V19 Cgg =0.93= Yy(c, Cc, 5) 


Fig 10—LiTaQ; stress-strain. 
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Fig. 11—Thermal stress comparison 2-growth LiTaQ; crystal. 


respectively. A plot of thermal stresses for the z-growth crystals by 
both methods is compared in Fig. 11. Agreement is very good. A 
comparison of thermal stresses by the long-cylinder method of both 
z-growth and y-growth is shown in Fig. 12 at @ = 45°. (Typically, 
6-stress variations are less than 5 percent.) The maximum stress 
recorded in absolute value is about 39,000 psi for both z-growth and 
y-growth. The z-growth (material axis 3 corresponds to x) would yield 
the same results as y-growth. This can be inferred from examination of 
the material properties of Fig. 10. 


VI. CONCLUSIONS 


Excellent comparisons were obtained for stresses and deflection in 
isotropic cylinders between the methods described in this paper and 
known theoretical solutions (Figs. 5 to 7). Stress comparisons again 
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Fig. 12—Thermal stress comparison z-growth and y-growth (at @6 = 45°) LiTaO; 
crystals—long-cylinder method. 


were excellent between the two methods described in this paper for 
z-growth LiTaO; crystals (isotropic in the plane of the circle) during 
the post-growth cooling stage (Fig. 11). Comparisons of thermal 
stresses between a y-(or x-) growth LiTaO; crystal (both anisotropic in 
the plane of the circle) and a z-growth crystal are presented (Fig. 12). 
Differences in the stress distributions were not great enough to favor 
elther growth direction; the more important consideration is to slow 
the cooling process down enough to keep radial thermal gradient to a 
minimum. 
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APPENDIX 
Stress-Strain Relationships—Long Cylinder 

The most general stress-strain relationship considered is the 4 X 4 
submatrix bounded by the dotted lines of eq. (18) as shown below, 
with no coupling between the stresses (2, cy, oz, Tzy) and the strains 
(Yue; Yen). 


Ox Giz C2, Ces 0 | 0 0 Ex 
Ty Cay. Cay Coe 0 10 0) Ie; 
Oz Cz; Ce Cx 0 | 0 O Ez 
Tzy( = 0 0 0 Ce. | 0 O Yeu’ (18) 
Ty: 0 0 0 0 V ye 
Tex 0 0 0 0 Vex 


After inclusion of thermal strains, eq. (18) can be restated as 


fox} = [Ce ta} ~ | faar}) (19) 


{oo} =a {oz Ty Fz Tau} 


{ €o} { €2 Ey €z Yeu}; 


| faar'| | fast foat foat 0}, 


and [C>] is the 4 X 4 submatrix of eq. (18). 

Now consider a rotation about z by the angle 6, and let {ea} 
= {o,0902T,} be the radial, tangential, longitudinal, and shear stress, 
respectively. The relationship between this rotated stress field and 
{oo} can be expressed as 


where 


I 


| 


{oa} = LB ]{oo}, (20) 
where 
cos76 sin?0 0 2 sin@ cosé 
sin24 cos26 QO —2 siné@ cosé 
Eee: 0 0 1 0 " 
—siné cos@é sinécos6 O cos?é — sin?é 
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Conjugate to eq. (20) is the following relationship between the strains 


{eo} = LB ]*{ea}, (20a) 
where 


| | { €a} ae { €, €9 €z Yro}. 
Kas. (19), (20), and (20a) can be combined as 
{oa} — [Ca ]{ ea} _ [B ][Co]{adT}, (21) 
where | . 
[C.J = [BIC JLB I. 
After multiplying eq. (21) by [C. ]—! and rearranging, the following 
result is obtained. 
{co} = [E]oa} + (CB) faa? (22) 
where 
(AJ = (C.J° = (B))-LO' ey. (22a) 


Tentatively, let e, = 0. This results in a longitudinal stress applied 
at the ends of the cylinder. From the third line of eq. (22), the longi- 
tudinal stress can be obtained as 


oz = 1/Bu( Ba Or ss E32 get Ess Tré => [oar), (23) 


where, from eq. (22a), 
E3; = H,,c0s?6 + Ey, sin’6, 
E32 = E,,s8in?0 + Hy, cos’6, 
Fi33 = Ezz, 

2 siné cosé (Ey, — Ezz), 

and E,,, E,., and E., are obtained from [ C5 |“. 


From eq. (23) we obtain 


a 
| 


{oa} = [D]{oo} — {00 1/E'33 faar of ; (24) 
where 
{oo} = {or 00 Tro} 
and 
1 0 0) 
0 ] 0 
pice 
LOI =) 8, /Bis By/Bis —Bu/Ba 
0 0 1 
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After substituting eq. (24) into eq. (22) and premultiplying by [D ]¢ 
(remember e, = 0), we obtain 


(<s} = [Es]los} + [Bs){ faa}, (25) 
where 
{eo} = { & €9 Yr0}, 
[Ao] = [DIL ILD] 
and 
[B,J = (DI(LB))7 | 
cos26 sin26 — E31/Es33 sin@ cosé 
= sin? cos? — E'32/E33 —sin@ cosé 
—2sin@cos@é 2sin0cos@ —LH34/E33 cos?@ — sin? 
Let [C] = LE, |~!. After inverting eq. (25) we obtain a desired result. 
jos} = [C1( {es} — [BsI| faa? ) (26) 


subject to the restriction ¢, = 0, which will be removed after an initial 
solution is obtained. 
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Mean-Squared-Error Equalization Using 
Manually Adjusted Equalizers 


By Y.-S. CHO 
(Manuscript received October 19, 1973) 


This paper describes an equalization procedure for systems using 
manually adjustable bump equalizers that 1s based on a mean-squared 
error criterion. We show that, 1n accordance with a Gauss-Seidel iteration 
process, gain adjustment always converges to the optimal value at which 
the minimum MSE of the equalized channel is obtained. Both zero 
forcing and MSE algorithms based on the Gauss-Seidel tteratton method 
are derived, and hardware implementation of these algorithms 1s dis- 
cussed. According to the error reduction analysis, an equalizer composed 
of orthogonal networks requires only one iteration to bring the equalizer 
to the opttmum state. For the bump equalizers used in the latest L5 Coaxial 
Carrier Transmission System whose Bode networks are semi-orthogonal, 
two to three iterations are shown to be sufficient to achieve the optimum 
gain settings in the mean-squared error sense. 


I. INTRODUCTION 


In this paper, a manual adjustment of bump equalizers is described 
which uses a mean-squared error (MSE) criterion. In the existing L4 
Coaxial Transmission System, the bump equalizers (realized with Bode 
equalizer networks‘) are used for line equalization and adjusted ac- 
cording to a zero forcing (ZF) algorithm.? This method results in an 
optimum equalization in the MSE sense, but only under certain very 
restrictive conditions with respect to the transmission response of the 
channel. The latest L5 Coaxial Transmission System, which provides 
up to 10,800 toll-grade long-haul message channels on a pair of 0.375- 
inch coaxial cables over 4000 miles, also includes bump equalizers for 
equalization of the 65-MHz bandwidth channel. The equalization 
method used in the L5 Coaxial Carrier Transmission System, however, 
minimizes the MSE of the channel deviation. It is found in practice 

847 


that the MSE algorithm has given a better equalization result than 
the ZF algorithm. | 

Since a coaxial cable system shows relatively stable channel charac- 
teristic, the bulk of the L5 line equalization is accomplished by 
manually adjustable equalizers. Normally, the time-varying channel 
deviation in a cable system is mostly the result of seasonal tempera- 
ture variation and aging of the components in the system. In such case, 
a usage of complex automatic equalizers in the system 1s not economi- 
cally desirable. 

In Ref. 3, two MSE methods are discussed for the equalizer adjust- 
ment. Both methods are based on the steepest descent algorithm 
and could be easily implemented in an automatic equalizer, but not in 
a manual equalizer. 

In this paper, an MSE algorithm based on the Gauss-Seidel itera- 
tion method is described for the gain adjustment of the manual bump 
equalizers. Under the specified assumptions, this method guarantees 
the convergence of the following iterative process. If a visual display of 
the gradient of the MSE with respect to each gain setting is available, 
adjust the first gain setting until the gradient becomes zero; next, 
adjust the second gain setting until the associated gradient becomes 
zero; similarly, adjust the third gain setting and all others up to the 
last one, thus completing one iteration. As the number of iterations 
increases, the residual error in the channel will be minimized in the 
MSE sense. 

While the Gauss-Seidel iteration method may seem quite compli- 
cated, it has several distinct advantages. First of all, the Gauss- 
Seidel iteration method requires only one gradient at a time, which 
can simplify the hardware involvement, particularly for manual 
equalizers. As is shown in Section III, the number of iterations needed 
to bring the equalizers to the optimum state is not large. When the 
equalizer is composed of orthogonal networks, a single iteration is 
sufficient. Since most of the equalizer networks used in transmission 
systems are orthogonal or semi-orthogonal in nature, the number of 
iterations will usually be small. For the bump equalizers, which consist 
of semi-orthogonal terms, two to three iterations are satisfactory for 
the optimal equalization according to the error analysis given in 
Section III. This result has been verified experimentally in the field. 


li. MSE ADJUSTING ALGORITHMS FOR BUMP EQUALIZERS 


In this section, several assumptions are made before the ZF and 
MSE algorithms are presented. 
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2.1 Characterization of coaxial channel 


The channel assumed in this section is represented by an infinite 
sin 2/x series on the frequency domain. (Since the transfer function of 
the Bode network is symmetric on the log f plane, where f is the 
natural frequency in hertz, the frequency used throughout is defined 
by w = log f.) 

Let M(w, #) represent the time-varying channel misalignment 
which is a real valued function of frequency in decibels. From the 
practical point of view, however, the channel can be assumed to be 
simply M(w) since the time variation is negligible during the equali- 
zation interval. Further, assume that the Fourier transform of the 
channel is limited in the time domain by a certain positive constant. 
(It should be noted that the Fourier transform of the channel does 
not result in an impulse response of the channel because the channel 
is measured in dB. However, there is an implicit dual relationship 
between the channel studied in this paper and that involving a time- 
domain equalizer, e.g., a transversal equalizer, in which a frequency 
band limitation of the channel is implied.) Hence, the channel can be 
characterized on the frequency domain by the following series: 


M(w) = ¥ Cc, Maple — wa)! (@B), (1) 


where C,, p, and w, are real numbers and 


Wnti — Wa = for alln = 0,1, -->. 


i 
2p 
Note that w, is equally spaced. 

Equation (1) can also be written in the following way. 


1 
M(w) = i; 2X, C',, cos[2rp(w — wn)x |dx 
0 n= 
} roe) 
= i X C', cos(2rpw,x) cos (2rpw) 
0 r= 
+ = C, sin(2rpw,2) sin(2rpwx) > dx 
1 
= [ {F(x) cos(2rpwx) + H(x) sin(2rpwz) } dz, (2) 
where 


F(z) = > C, cos(2rpw,x) and A(z) = 3 C, sin(2rpw,t). 
=0 =O 


nm nr 
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Since 0 S x S 1, eq. (2) implies that the shortest frequency domain 
ripple period found in the channel M(w) is 1/p. 


2.2 Representation of bump equalizers 


The bump equalizer considered in this paper is a linear combination 
of adjustable-loss Bode networks. The input-output transfer function 
of an equalizer composed of N Bode networks can be represented by 


EQL(w) = ¥ geBe(w) (AB), (3) 


where N is the number of networks and g; and B, represent the gain 
and response respectively of the kth Bode network. 

A typical Bode network is shown in Fig. la, where the loss is con- 
trolled by the resistor R. The transfer function, B;(w), can be analyti- 
cally derived and, with a suitable flat gain amplifier, it can be expressed 
by the following equation: 

= [£01 + By) + D, |? — D, 


B,(w) = ~ (G+ Ey? + Die (dB), (4) 


where 
Ror 

kK, == 

ae 

D,z _ (w/wz) Ai 

* [(w/wi)? — 17?’ 
_ Cx 

H,. = ie 

and 


1 
a ee 
ee eS QavlinC;, 


Since eq. (4) shows B;(w) to be a quite complicated function of w, 
one of the following assumptions is used while analyzing the equalizer 
in detail. 


Assumption 1: Let B;(w) be approximated by 


sin[a(w — wz)/Aw ] 


= ese Tae (dB). (5) 


. vis 
sine| Aw — Wz) | = 
Since there are N Bode networks in the equalizer, which for the 
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Fig. 1—Adjustable Bode network. 


present analysis are spaced equally on the w-axis with interval Aw 
(see Fig. 2), then the transfer function of equalizer can be expressed by 


N’ T 
EQL(w) = Py rE sine| Aw! — wz) | (dB). (6) 
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Fig. 2—Manual equalization of L5 coaxial system. 


To permit a comparison between B;,(w), as represented by eqs. (4) 
and (5), the two equations are plotted in Fig. 1b. The maximum dif- 
ferences between the two best matched curves are 0.165 and 0.183 
dB when |w — w;| S Aw and |w — w;| > Aw, respectively. 


Assumption 2: Let B,(w) be approximated by 
_ sinLr(w — wx)/Aw] 
a(w — wz)/Aw 


— cos[r(w — wz) /Aw | (7) 
1 — 4[(w — w;z)/Aw J? 


cosinc| —-(w — wz) 
Aw ; 


Under the same conditions listed in Assumption 1, the transfer func- 
tion of equalizer can be expressed by 


EQL(w) = = Tk cosine| rence — wz) | (dB). (8) 


Expression (7) is also plotted in Fig. 1b, and it can be seen that cosine 
[r(w — wx)/Aw] approximates quite well the actual transfer function 
of the Bode network as expressed by eq. (4). The maximum differ- 
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ence between the two best-matched curves is 0.0327 dB when |w — w; 
< Aw and 0.0404 dB for |w — wz] > Aw. 


2.3 Mean-squared error algorithms (Gauss-Seidel iteration method) 


The definition of optimal equalization used in this paper is the mini- 
mization of the MSE of the equalized channel as a result of adjusting 
the gain parameters, gx. 

On a decibel scale, the equalized error will be 


N 
E(w) = Pp» gi.Bi(w) — M(w). (9) 
Then the MSE can be represented in the frequency domain by 
MSE = / ” 1 E(w) | 2d. (10) 


Theorem 1: If the equalizer described by eq. (8) is composed of linearly 
independent networks, then there exists a unique set of gi’s which nulls 
all the gradients Gy; 2.€., 


_ aMSE 


Gi aa 


= 0 forallk =1,2,---,N, (11) 
where MSE is defined in eq. (10), and the corresponding set of g:’s 
results in mintmum MSE. 


The proof is given in Appendix A. 

The bump equalizers considered in this paper belong to the class 
of the equalizers defined in Theorem 1. 

As derived in eq. (25) of Appendix A, the gradient vector is given by 


G = Bg — M, (12) 


where B, g, and M are the system matrix, gain vector, and correlation 
vector, respectively, and are defined as follows. Defining an inner 
product 


(A, B) = is A(w)B(w)dw, 
G= LG, Ga, my Gy I’, 
where 7' indicates the transpose, 


= Lo, G2, °"°, gv |’, 
M = 2[ (Bi(w), M(w)), (B.(w), M(w)), ees (By(w), M(w)) J’, 
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and 


(Bi, B;), (By, Ba»), an (Bi, By) 
B = 2 (Ba, By), (Bo, Bz»), “——™ (Bo, By) 
(By, By), (By, Bz), or (By, By) 


Now the gain vector is 
g = B-(G + M), (13) 


provided that B— exists. The optimum gain, g*, which results in the 
minimum MSE is obtained by solving eq. (13) with G = 0. The 
MSE algorithm given in Ref. 3 solves (13) with G = 0 by the steepest 
descent method, which can be readily implemented in an automatic 
equalizer control circuit. For manually adjusted equalizers, however, 
the steepest descent algorithm cannot be easily implemented because 
the algorithm requires simultaneous adjustment of all the gain settings. 
A manual equalizer adjustment algorithm should have the following 
properties: 


(2) The several gain settings can be adjusted one at a time, within 
a specified sequence. 
(22) Repeating step (2), g approaches g*. 


The converging rate of the initial g to the final g* depends on the type 
of algorithm used and the system matrix B. For the bump equalizer, 
this question is discussed in Section 3.2. 


Theorem 2 (Gauss-Seidel iteration algorithm): If the system matrix 
B in (12) has dominant diagonal elements such that 


|(By Bs)| > Zt |(Bs, Bs)| (14) 


for all k=1, 2, ---, N, 


where >’ indicates the summation of all terms excluding the case 7 = k, 
then every g converges to the optimum gain g* = BM by the following 
iteration process: 


Iteration 1: Let gzcy indicate the kth gain at the ith iteration; thus, the 
initial gain settings are gio), 9200), °° *, JN). Adjust gio) until ats corre- 
sponding gradient G, = 0 and designate the resultant gain gia). The gain 
settings are then gia), 920), *°*, Jno). Adjust goin) until the gradient 
G2 = 0, resulting in the gain settings gia), Gea), 930), °*°) JNO): 
Repeating the operation for each setting results an gic), Joa), 93a), °°; 
gn1), and completes the first iteration. 
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Iteration 2: Adjust gia) until Gi = 0, resulting in the gain settings gic), 
Jay, 93), ***, YN). Obtain gaia), 93(2), +++, Jn (2) by similar operation, 
completing the second iteration. Similarly, tterations 3, 4, --+, n can be 
carried out as requtred. 

The proof is given in Appendix B. 

If equalizer networks satisfy the inequality (14), they are called, 
in this paper, semi-orthogonal terms. The bump equalizers defined in 
this paper satisfy the inequality (14), and hence Theorem 2 can be 
used as a manual-equalizer adjusting algorithm. To implement this 
algorithm, a visual display of each gradient is required before the cor- 
responding gain is adjusted. The following two theorems provide 
simple ways of determining the gradient. 


Theorem 3 (ZF algorithm): Let the channel be represented by eq. (1) and 
the equalizer satisfy assumption 1. If the interval Aw between two adjacent 
Bode networks is no greater than half the shortest ripple period ( = 1/p) 
found in the channel, 1.e., 

(15) 


Aw 2p 


IA 


then the optimum gain setting is obtained by repeating the Gauss-Seidel 
iteration process defined in Theorem 2 with the gradient given by 


G, = 2H (wz), (16) 


where k = 1, 2, ---, N and E(w,) is the frequency domain error value 
measured at frequency wz, the center frequency of the kth Bode network. 

The proof is given in Appendix C. 

Thus, if signals are transmitted at a set of frequencies equal to the 
center frequencies, w;, of the Bode networks, and if the errors are 
measured at these frequencies at the receiving station, the gradients, 
G;,, can be obtained directly. Then each gain, g;, would be adjusted 
until its gradient, G., reduced to zero. This is the well-known ‘‘zero- 
forcing’ technique used in Ref. 2; it also achieves the optimum equal- 
ized channel in the MSE sense, if the stipulated assumptions apply. 


Theorem 4 (MSE algorithm): Let the bump equalizer in this case satisfy 
assumption 2 and assume that the interval, Aw, between adjacent Bode 
networks 1s no greater than the shortest frequency domain ripple period in 
the channel, 2.€., 


Aw 


IIA 
<S lee 


(17) 
Then the optimum gain setting 1s obtained by repeating the Gauss-Seidel 
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iteration process with the gradient in this instance given by 
A 
k= 1, 2, neg dN, 


where E(w,;) ts the frequency domain error at the center frequency of 
the kth Bode network, and E(w; — (Aw/2) |] and H[w, + (A/2) | are the 
frequency domain errors measured at lower and upper frequencies midway 
between adjacent Bode networks. Equation (18) 1s derived in Ref. 8. 


Proof: In this case, 


and 
N 


do’ |(Bx, By)| = 0.25 


j=1 


for all k, thus satisfying the inequality (14). Hence, the Gauss-Seidel 
iteration process converges to the optimum gain settings. 

To implement the MSE algorithm, a measure of the error at 2N — 1 
points in the frequency domain is required (see Fig. 2). In practice, 
the MSE technique results in better equalization than that obtained 
by the ZF method. Note that assumption 2 for the MSE algorithm 
approximates the actual equalizer more precisely than assumption 1 
does. Moreover, inequality (15) for the ZF algorithm derived in this 
section is a conservative assumption. The channel ripple period allowed 
by the MSE algorithm can be half the period assumed by the ZF 
algorithm. 


lll. CONTROL OF MANUAL BUMP EQUALIZERS 


In this section, the Gauss-Seidel iteration process derived in the 
previous section is applied to the manual equalizer for optimum gain 
control. The number of iterations required to obtain acceptable gain 
settings is reflected in inequality (14). The larger the diagonal com- 
ponents [left-hand side of (14) ] compared to the off-diagonal com- 
ponents [right-hand side of (14) ], the fewer the iterations needed. The 
rate of convergence of the iteration is described in Section 3.2 based on 
an error-reduction analysis. It is shown that one iteration is sufficient 
to obtain the optimum gain settings for the ZF Gauss-Seidel iteration 
algorithm derived in Theorem 3. When, in the more general case, the 
channel is initially equalized by the ZF algorithm, one or two more 
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iterations are usually sufficient to achieve a practically optimum 
equalized channel for the MSE algorithm. | 


3.1 Hardware realization of Gauss-Seidel iteration process 


For the L5 line equalization, an equalizer adjustment system has 
been developed for the adjustment of bump equalizers by the Gauss- 
Seidel iteration process. It 1s composed of a precision transmission 
measuring set, 90G oscillator—90H detector—digital control unit 
(Ref. 4), and a hardwired, special-purpose computer which contains 
a programmed memory and an arithmetic unit called an equalizer 
adjustment unit (HAU) (see Fig. 3). Referring to Fig. 2, we assume 
that the equalizers in the receiving station are to be adjusted by the 
MSE algorithm. By selecting the particular Bode network to be 
adjusted, the EAU in the transmitting station causes the 90G oscillator 
to generate sequentially an appropriate set of three frequencies (fi1, 
fio, and fz3). The KAU in the receiving station measures the channel 
error at the same three frequencies and computes the gradient, Gz, 





Fig. 3—Equalizer adjustment unit. 
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by the following relationship : 
G, = BirE (fir) + Broki( fre) + Byski( fis). 


Normally, By, = Bis = 3 and By. = 1. The resultant is displayed in 
a digital readout of the EAU and the operator subsequently adjusts 
g, until G;, = 0. Then the next Bode network is selected and the trans- 
mitting EAU causes the generation of the required set of three fre- 
quencies, and the receiving EAU processes the received error signals 
and displays the calculated gradient. Again, the corresponding gain 
adjustment is made. When the ZF algorithm is selected, the gradient 
displayed is simply the error at the center frequency of the Bode 
network. Hence, as far as the operator is concerned, the adjustment 
procedure for the ZF and MSE algorithms is identical. 

In practice, it is found that the equalizer should be initially adjusted 
by the ZF algorithm to bring the system near the optimum state. In 
this way, any large initial gain deviations from the optimum value are 
quickly reduced to within about +0.5 dB. Then the MSE algorithm is 
used for the ‘‘fine tuning’ of the gain adjustments. Usually, one or 
two iterations with the MSE algorithm will be sufficient for the equal- 
izer to reach the optimum MSE state when starting from the ZF 
state. According to inequality (23) derived in the following section, and 
given initial gain settings within +0.5 dB of the optimum value, the 
gain settings after two iterations are within +0.036 dB of the ideal 
values in the worst case. The actual deviation from optimum in most 
cases will be smaller than 0.036 dB and in any case will be less than 
the inherent accuracy limitations of the 90-type transmission measuring 
sets to be used. 


3.2 Error analysis 


The Gauss-Seidel iteration provides fast convergence of initial gain 
settings to the optimum values for the bump equalizers. As derived in 
Appendix B, the Gauss-Seidel iteration can be expressed by 


Sis1 = — L'Uga, + LM, (19) 


where L is the lower triangular portion of the system matrix B in- 
cluding the diagonal and U is the upper triangular portion of B not 
including the diagonal. 

Defining an error vector after the zth iteration to be 


Qu) = £* — Za), (20) 
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where g* is the optimum value, then 


C41) =B* — B41) 
= g* + L-'Ugi) — L-'M 
= L~Ue(:). (21) 


To derive (21), the equality, L-'M = g* + L—!Ug*, was used. Hence, 
the magnitude of the eigenvalues of L~'U determines the speed of 
convergence. 

If the equalizer belongs to the class defined in assumption 1, Bis 
a positive definite diagonal matrix and L~'U is a null matrix. Hence, 
Qu41 = Oforallz = 0,1, 2,---. Inother words, we obtain the optimum 
gain settings by the ZF Gauss-Seidel iteration algorithm in one iter- 
ation. The same result can be obtained from eq. (19), since L~!U is a 
null matrix and L7! = B-, 

If the equalizer satisfies assumption 2, the MSE Gauss-Seidel 
iteration algorithm can be used. Now the system matrix is 


641, by, bi3, sey bin 
B = ») bat, boo, bos, aaa bey 
bw1, be, bnx, +++, Onn 
where 
b;; = 0.75 if ~= 7 
and 
b;; = 0 Meg ee 
for all 2, 7 = 1, 2, ---, N. Splitting B into two parts, L and U, which 
are defined above, and performing some algebra, 
0, 67, 0, O, ---, O 
0, —6-%, 67}, 0, ---, O 
L0 = 4|0, 67%, -67%, 671, ---, O 
0, (—1) 4416-4, (—1)%6-@-), ae 6-3, — 6-2 


Hence, one can calculate the new error vector by eq. (21). After one 
iteration, the upper bound on the maximum residual error becomes 


| €x (1) | max = 3(67? = 6-2 ++ 6-8 + nae -) | €¢0) | max 
= 0.26667 | €5(0) | max (22) 


for all j, k = 1, 2, ---, N. Similarly, after the second iteration, 
| €x (2) | max <= 0.26667 | C5 (1) | fase 0.07111 | @z(0) | max (23) 
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for all 7, j, k = 1, 2, ---, N. The equality in eqs. (22) and (23) will 
be obtained if and only if 


—~ — —_— 


C10) = = 20)" = C30) = — 2410) Se Crt: wines 


This in general will not be the case, and the maximum setting error 
after the second iteration usually will be less than 0.07111 times the 
maximum setting error prior to the first iteration. Consequently, if 
the gain settings are within 0.5 dB of optimum at the start (which a 
single ZF iteration will establish), the deviation from optimum settings 
in the MSE sense is within a few hundredths of a decibel after two 
additional iterations. 


IV. CONCLUSION 


This paper shows that a manual equalization process which can be 
described by a Gauss-Seidel iteration method provides optimal con- 
trol for bump equalizers in the MSE sense. Compared to the steepest 
descent method discussed in Ref. 3, the Gauss-Seidel iteration method 
can be more economically implemented for manual equalization such as 
in the L5 system. The Gauss-Seidel iteration process requires knowl- 
edge of the gradient of the MSE with respect to the gain setting for 
each Bode network to be adjusted. The ZF algorithm derived in this 
paper requires just one Gauss-Seidel iteration, but in practice the ZF 
equalized channel is not optimum and can be further improved by the 
MSE algorithm. This is because the gradient obtained by the MSE 
algorithm is more accurate than the one obtained by the ZF algorithm 
for realizable equalizer shapes. It should be noted that three tones are 
required to determine the gradient of the MSE with respect to the 
gain setting of each Bode network for the MSE algorithm, while only 
one tone is used to obtain the gradient for the ZF algorithm. The 
number of iterations that are necessary to bring the equalizer to the 
optimum state depends on how close the initial settings are to optimum. 
When the channel is initially ZF-equalized, only one or two more 
iterations are needed to optimize the channel with the MSE algo- 
rithm. This result agrees completely with experiments conducted in 
the Ld field trial. 


APPENDIX A 
Proof of Theorem 1 


Substituting eqs. (9) and (10) into eq. (11), we obtain 


G,=2 1] B,(w) a 9;B,w) — M(w)\ dw. (24) 


=O 
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Defining an inner product 
(A,B) = [~ A(w)B(w)du, 
pote nea: 
where 7' indicates the transpose, 


g = (91, 92, ++, gn]? 
M = 2[(Bi(w), M(w)), (B2(w), M(w)), Ts (By(w), M(w)) |? 


and 


(Bi, Bi), (Bi, B,), eg (Bi, By) 
B = 2 (Bs, B,), (Bo, Bo) , mete ey (Ba, By) 
(By, B;), (By, Be), aaa (By, By) 


a simultaneous equation of the type of eq. (24) for all k from 1 to N 
can be written as 

G=Bge—M 
or 

Bg = G+ M. (25) 


Since eq. (25) is a nonhomogeneous system of N equations and G + M 
is a vector with N real-numbered components in the case considered, 
for the given G, the unknown g is uniquely obtained by 


g=B(G+ M), (26) 


provided that B is a nonsingular matrix. 
However, if B were a singular matrix, then a linear combination of 
the columns could be made zero, i.e., 


N 
2 hi Bx, B;) — 0 


for each j = 1, 2, ---, N where h, are real nonzero numbers. 
In addition, the following relationship could also hold: 


N’ N 
Py Ayhy( By, Bi) + p> hehi( By, Be) + ++ 


N 
+ Py hyh,(B;, By) = 0. (27) 
But (27) can also be written as 


fe ae heBa(w)} du =i6, (28) 
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Equation (28) contradicts the assumption that B,(w)’s are linearly 
independent. Hence, B is indeed a nonsingular matrix and there exists 
but one set of g,’s for which G = 0 in eq. (26). It is yet necessary to 
prove that this stationary point is the global minimum of the MSE 
defined in eq. (10), which is established if the following relationship 
can be proved: 


[| arc) = aX oiBatw) ~ 8 Xo Bilw)} aw 
/ ro f N ; )2 
Ce fi Mele) — 2 geBalwo) | dw 


00 N 2 
+e | {aw _¥$ rB(w)| dw, (29) 
= k=l 


where 
x * x ek 
ps, g,.B,(w) and Py g; By (w) 


indicate distinct equalizer settings, a+ @= 1 and a@ > 0, then 
MSE is a strict convex function of gain settings g;’s and has a global 
minimum. 

Subtracting the left-hand side from the right-hand side of inequality 
(29), we obtain the following: 


— 2a a i p> GiB (w) | + = Bu) || dw. (30) 


Since B,(w)’s are linearly independent and at least one equalizer 
setting, 


N N 
uu grBi(w) or 2 gx Bi(w), 


is not zero, (30) is negative. Hence, inequality (29) is correct, and the 
proof of Theorem 1 is complete. . 


APPENDIX B 
Proof of Theorem 2 


_ The gradient of MSE with respect to the gain settings is represented 
by the following equation: 


G = Bg — M, (31) 
where B, g, and M are defined in (24). Splitting the B matrix as follows 
B=L+U, (32) 
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where L is the lower triangular portion of B including the diagonal and 
U is the remainder of B, 


G = Lg + Ug — M. (33) 


According to the iterative procedure described in the theorem, g (41) 
is obtained from g,;) by setting G = 0. With the aid of eq. (33), this 
procedure can be expressed by 
= Lgcsi) + Uga, — M 
or 
Susy = — LUga) + LM. (34) 


By successive calculation, eq. (84) can be modified by 
Gory = (-L IU] g¢0) + 2 [-L-U LM, (35) 


where g (0) 18 the initial value. 
If 
[-L-Uu]'—[0] as t> o@, 


Se P-L-UyL-! > [L + Up! = Bo. 
k=0 


Hence, eq. (35) becomes 
Zu+1 = O+ BM, 


which is the desired result. | 

Hence, the theorem is proved if L~!U is a convergent matrix, i.e., 
eigenvalues of the matrix L~!U are all less than one in absolute value. 
However, if condition (14) is satisfied, L~'U is a convergent matrix 
and [—L~U ]' > [0] asi — ~ (see Theorems 3.3 and 3.4 in Ref. 5). 


Note: The iteration process defined by (384) is known as the Gauss- 
Seidel or, simply, the Seidel iteration. 


APPENDIX C 
Proof of Theorem 3 
When assumption 1 is satisfied, 
(B,, By) = 1, k= J 
=0, k¥] 
for all 7 and k. 
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Hence, Theorem 2 is satisfied and we have to prove now 


CG = 2E (wx). (36) 
From eq. (24), 


00 N 
G, = 2 [ Ba(w)| p> g;B;(w) — M(w) jaw 
=2 /. sine| Aj (w ee) pe g;B;(w)dw 
eS, iZ sine| x, (w — wy) | (w)dw (37) 
Since * 


[ sine| rene — uw) | sine| Ate — w) fav =0 if k#j 


and 
1 if k=4j, 


the first integration in (37) is simply 2g9,. 
Substituting M(w) of (2) into (87), the second integration of (37) 
becomes 


2 [ sine exc — wi) | f { F(x) cos(2rpw2) 
: + H(z) sin (2rpwx)}dxdw 
= [~ sine | = wi) | [ (4@) cosl2xp(w — w.)2} 


+ h(x) sin [2rp(w — wz)x|}drdw, (88) 


where 

f(x) = F(x) cos (2rpw,x) + A(x) sin (2rpw;c) 
and 

h(x) = H(z) cos (Qrpw,x) — F(x) sin (2rpw,2). 
Since 


1 co 
i | h(x) sine | Tu — wi) | sin [2rp(w — w,;)x |dwdx = 0, 
0 —o 
Eq. (38) becomes 


2 [ sine | Zw — wi) | ia { f(x) cos [2rp(w — wz)x]}drdw. (39) 


) 


Replacing w = u + w, and changing the ‘‘cos” into ‘exponential’ 
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form, (89) becomes 


J ° sinc (=) [ { f(z) Lexp(22rpux) + exp (—722rpux) |} drdu 


. Os A T 
= [ f(x) be sine (4) {exp(12rpuz) 
+ exp(—72rpux)}dudz, (40) 
where 722 = — 1. Since OS x S 1 and 2p S 1/Aw by assumption, 


integration of (40) is simply 


2. i: ade. 


Note that the inner integration in (40) is the Fourier transformation of 
the sine function. 
Combining the results, G;, in (87) becomes 


a ee: [ coe (41) 


However, 
1 
M(w,) = i! sie wae BOLGy- =a. 


Hence, (41) becomes 
G, = 2LEQL(w,) — M(wz) ] 
= 2H (wz). 


This proves the theorem. 
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Mathematical Analysis of an 
Adaptive Quantizer 
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This paper presents a mathematical analysis of an adaptive quantizer, 
a pulse code modulator, which 1s used for coding speech and other continu- 
ous signals with a large dynamic range into digital form. The device is a 
two-bit quantizer in which the step size is modified at every sampling instant 
with the object of adapting the range of the device to the intensity level of 
the signal. In the adaptation algorithm analyzed in the paper, the encoded 
information of the previous sampling instant is used either to increase 
or to decrease the step size by fixed, but not necessarily equal, proportions. 

Initially, the stochastic stability of the device is established by construct- 
ing a stochastic Liapunov function. Various basic identities and bounds 
on aspects of the behavior of the device are obtained. The qualitative 
results obtained indicate the nature of the trade-offs between the quality 
of the steady state and the transient performance of the device. Also, 
formulas are developed for the purpose of evaluating the mean time 
required for the step size to adapt from arbitrary initial conditions to 
certain optimal values. 


I. INTRODUCTION 


A mathematical analysis of an adaptive quantizer is presented in this 
paper. The coding thresholds of the device, also referred to as the step 
sizes, are not fixed but adapt according to a particular alogrithm. The 
object of the algorithm is to modify the threshold to larger or smaller 
levels, depending on whether the signal intensity level is high or low, in 
a manner that allows a decoder at the receiving end to effectively re- 
construct the continuous signal. The basic two-bit quantizer, 1.e., quanti- 
zers with four output levels with codes 01, 00, 10, and 11, is character- 
ized by a particular function of the following form at each sampling 
instant. 
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Input refers to the nth sample of the continuous signal, x(n), 
n = 0,1, 2, ---; output refers to the coded signal to be transmitted at 
that instant; and A is the step size. In adaptive quantizers of the type 
to be investigated here, the step size is variable and the step size at the 
nth sampling instant is denoted by A(n). The step size uniquely defines 
the entire function in the manner indicated by Fig. 1; hence, the com- 
plete adaptive quantizer is associated with a sequence of functions. 
The adaptive quantizers that are the subject of this paper are basically 
characterized by the following adaptation algorithm 


A(n +1) = M,A(n) if |a(n)| S A(~m) (la) 
= M.A(n) if |x(n)| > A(n), (1b) 


where M, and M,, called multiplier coefficients, are fixed constants 
satisfying* 0 < M, <1 < Mz. Variations on (1) are considered in 
the main text, although the discussion in the introductory section is in 
terms of (1). Results on adaptive quantizers with output levels more 
numerous than 4 will be considered in a future publication. 

The adaptation algorithm in (1) is due to Cummiskey, Flanagan, 
and Jayant.!? In Ref. 1 Jayant presents the results of extensive com- 
puter simulations undertaken to determine the multiplier coefficients 
which maximize various performance functionals. A class of random 
inputs {x(n)} that is considered is obtained by passing a discrete, 
white, Gaussian process through a filter with a single pole. In Ref. 2, 
Cummiskey, Jayant, and Flanagan consider a differential PCM coder 
in which the adaptive quantizer is used together with a fixed first- 
order predictor in the feedback loop. Their work has its direct ante- 
cedents in the various schemes?:*> for adapting step sizes in delta- 
modulators, a one-bit quantizer, and in the work of Wilkinson.® Wilkin- 
son’s paper on a two-bit adaptive quantizer, largely concerned with 
hardware implementation, is particularly interesting. In his scheme, the 
step size is controlled by a moving fraction obtained by keeping a tally 
of the number of times the input falls in the lower slot of the quantizer. 
Goodman and Gersho’ have independently looked at the adaptive 
quantizer from a theoretical standpoint and their work complements 
the work described here. 

In this paper we make a number of simplifying assumptions about 
the input sequence {x(n)}, the most restrictive being the assumption 


“Since the absolute value of the input in Fig. 1 is partitioned into [0, A] and 
(A, » ], we shall loosely refer to the event leading to (1a) as “‘the input falling in the 
lower slot.”’ 
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OUTPUT 





Fig. 1—The quantizer function. 


that it is a sequence of independent random variables. However, we 
have obtained for the idealized model precise results which indicate 
rather fully the trade-offs involved in the choice of the multiplier 
coefficients. Also, we have developed formulas for efficiently computing 
functionals as aids in the design problem. We believe that the broad 
qualitative features of the device that are found to hold in this model 
carry over for more realistic input processes. It is hoped too that the 
techniques developed here will provide a point of reference for future 
work. 

The mathematical analysis, for the main part, is of a random walk 
on the integers, whose complexity is due to the dependence of the state 
transition probabilities on the states. The structure of the random walk 
which is exploited here is rather general, and for this reason the model 
is of independent interest; to our knowledge, the main mathematical 
results have not appeared in the literature on random walks. 

The organization of the paper is as follows. In Section 1.1 we con- 
tinue the discussion on the adaptation algorithm in the context of a 
particular idealized model of the sequence {x(n)}, and we discuss some 
of the results to be derived later and what is already known about 
optimal quantization in the nonadaptive framework. In Sections 1.2 
and 1.3 we give the basic equations of the process arising from (1), 
and certain modifications of it, when the input sequence {z(n)} is 
independent and identically distributed. In Section II the stochastic 
stability of the device is established under general conditions. The 
existence and uniqueness of the stationary distribution of the step 
size is proved by constructing a stochastic Liapunov function for the 
random process. Section III examines in detail the stationary step 
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size distribution. In Section 3.2 we prove an identity which explicitly 
gives the stationary probability of the input falling in the lower slot 
of the quantizer, i.c., Pr. [|x(n)| S A(n)]. In Section 3.3 sharp 
bounds are obtained on the stationary probabilities. It is shown that 
for almost all values of the multiplier coefficients there exists a natural 
center of the distribution and that the stationary probabilities fall off 
at least geometrically with increasing distances from the natural 
center. In Section 3.5 results are obtained on a particular limiting 
behavior, namely, the effect of the stationary distribution of making 
both multiplier coefficients close to unity. Section IV is devoted to the 
transient response of the device. In Section 4.1 we develop formulas 
for the efficient computation of the time required for the step size to 
adapt from an arbitrary initial value to the desired step size. Section 
4.2 by giving an explicit bound on this time provides some insight into 
the dependence of the adaptive time on the choice of the multiplier 
coefficients. Finally, we report some computational results. 


1.1 Background 


In an idealized model for the samples, x(n), of the continuous signal 
process, assume that {x(n)} is a sequence of independent random vari- 
ables with zero mean. Assume further that the distribution of x(n) 
for every n is an element of the same equivalence class of distributions 
in which the distributions are equivalent to within a scaling operation. 
The scaling or intensity level changes slowly with n. For instance, the 
equivalence class of distributions may be the family of Gaussian 
distributions and only the variance, indicating the intensity level, 
changes with n. 

It is necessary to recall at this stage some known facts concerning 
the design of quantizers in the nonadaptive framework’ where {x(n)} 
is a sequence of independent, identically distributed random variables 
and the step size is fixed. Suppose that E[{y(n) — «(n)}?] measures 
the performance of the quantizer where y(n) is the nth output of the 
device.* The step size which minimizes this functional, A, is in principle 
easy to establish, and A is uniquely characterized by the probability 
of the input falling in the lower slot, ie., Pr[|xz(n)| < A]. Another 
observation that is equally easy to verify is that the optimal step size 
has the property that if the distribution of {2(n)} is scaled, then the 
optimal step size is obtained by an identical scaling of the previous 
optimal step size. A convenient way of stating this observation is: a 


“It is not essential that the performance functional be of that form. 


870 = THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1974 


property of the optimal step size that is invariant to scaling of the 
distribution of {x(n)} is the probability that the absolute value of the 
input x(n) does not exceed the optimal step size. For instance, when the 
distribution is Gaussian it is known that this probability is close to 
0.68.8 

An intermediate step in proceeding from the nonadaptive case to the 
more general model described prior to it, in which the identically 
distributed condition does not hold, is provided by the following model. 
Assume that the sequence {x(m) } is indeed independent and identically 
distributed, and that the equivalence class of distributions to which 
the particular distribution belongs is known. However, the scaling 
parameter is unknown. It is relatively straightforward to state the 
requirements on a well-behaved algorithm operating in this simple 
framework, and, if these requirements are always satisfied, then it is 
possible to conclude that the device will operate satisfactorily for the 
more general model. The requirements are: (2) for arbitrary initial 
step size guesses, the step size rapidly converges to the optimal step 
size, and (72) it is thereafter localized in a small neighborhood of that 
point. This paper separately analyzes the two requirements in the 
simple framework just described. Considerations related to (2) and 
(12) are lumped respectively under the terms ‘‘transient response” and 
“‘steady-state response,”’ since the latter property is effectively investi- 
gated in terms of the stationary distribution of the step size, assuming 
one exists. A good reason for the division is that they lead, in some 
ways, to quite opposite requirements for the multiplier coefficients. 

Consider, in the light of what is known about optimal quantization 
in the nonadaptive framework, what is required for the localization 
property, requirement (77), to hold. When the stationary distribution 
has both of the following properties, it is possible to establish an effec- 
tive correspondence and infer that (zz) holds: (a) the stationary proba- 
bility of the step size falling in the lower slot, i.e., Prs [|x(m)| S A] 
equals the known value associated with the particular family of dis- 
tributions; and (b) the mass of the stationary distribution is concen- 
trated in the small neighborhood of a point. In Section III we show that 
by appropriate choice of the multiplier coefficients it is possible to 
achieve both requirements. 


1.2 Basic assumptions and equations 


We consider only quantizers with multiplier coefficients having the 


following structure: 
M,=y7-" and M,.= y’, (2) 
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where y is some real number greater than 1 and k and IJ are positive 
integers. We shall further make & and | relatively prime, i.e., their 
greatest common factor is 1. If, as we shall assume, the initial step size 
is of the form y*, with ¢ an integer, then the step size is always of that 
form and the space of possible step sizes forms a lattice.* 

There is a step size with, as we shall see, certain claims to being the 
central step size for a particular distribution of {a(n)} and choice of 
parameters k and 1; this step size is used as a reference point. There 
exists an integer 7 such thatt 

PrCio@|Sv1< ps5 SPr[ie@l| S72 @) 
We denote y' by C and refer to it as the central step size; all step sizes 
are considered to be of the form Cy‘, 7 = 0, £1, +2, --: 

Obviously, it is more convenient to work with the log transform of 
the step size, so let 


w(n) & log, A(n) — log, C. (4) 


From the original algorithm we have 


a(n +1) = w(n)—k if [a(n)| S Cye™ 
=a(n) +l if [a(n)| > Cye™. (5) 
We have in (5) a Markov chain with states 0, +1, +2, ---. The state 


transition probabilities are obtained from the distribution of x(n): 
for all integers 2 let 
b; = Pr[|z(n)| S Cy‘] (6) 
and 
ie = 1 be 


The ‘‘b” is a mnemonic for backward probabilities since it is associated 
with a transition backwards from the generic state 2 to (¢ — k). The 
diagram in Fig. 2 represents the Markov chain. Denoting by p;(n) 
the probability that w(n) = 7, we have 


(7) 





Although the transition probabilities depend on the distribution of 
x(n), the two following properties of the sequence {b;}, on which we 


*D. J. Goodman suggested the above euuete on the multiplier coefficients with 
the object of obtaining a discrete Markov proces 
We are tacitly assuming that Pr [|z(n)| = 0) < <l/(k+) —e,e>0. 
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itd +1 





b; bi +1 


Fig. 2—The Markov chain. 


base our results, hold irrespective of the distribution: 


056; < bi S 1. for all i, (8) 


and 


bg eo es Sh (9) 


That the strict inequality in (8) holds for all 2 is a mild restriction on 
the distribution of x(n); however, certain straightforward modifications 
may be made to obtain corresponding results when the strict inequality 
does not hold for all 7. 

The property of the 0 state to which we alluded earlier may be loosely 
stated, thus: there is a net drift to the left (right) from states to the 
right (left) of the 0 state. Formally, 











Elw(n+1) |o(n) =1]—7= — (kK +) Cees | <0 if 7>0 
gs (10) 





>0 if <0. 





The above super- and submartingale properties are the basis for the 
existence of a stochastic Liapunov function (Section 2.2) and the 
bound obtained in Section 4.2. 


Remarks: The random walk in (5) with k = 1 = 1 is also the modei 
for the delta-modulator subject to random, independent, identically 
distributed inputs. The stationary behavior of the model was treated 
in an elegant paper by Fine.® Gersho” has established the stochastic 
stability of the delta-modulator for a larger class of input processes. 
Some of our results, particularly those in Section IV on transient re- 
sponse, appear to be new and of some interest in this context. 


1.3 The saturating adaptive quantizer 


For the algorithm in (1) and, say, Gaussian distributions of the 
input, there is a small, positive probability of the step size exceeding 
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any large prespecified level. A model which reflects more accurately 
the practical algorithm for adapting the step size is one which does not 
allow the step size to become unbounded. One way of implementing 
this is to make the step size saturate at some suitably large level, 
ie., if A(n) < |a(n)|, then 


A(n + 1) = min[M,A(n), L]; L> 0; (11) 
i.e., in the log transformed variables, 
w(n +1) = min[w(n) +10]; L>O. (12) 


The model of this device, which we shall refer to as the saturating 
adaptive quantizer, is useful not only for the reasons given but also 
on theoretical grounds since the results obtained for the saturating 
adaptive quantizer yield, in the limit as L > ©, corresponding results 
for the adaptive quantizer. We carry both models with us throughout 
the paper and at least indicate along the way the main correspondences. 

Tor similar reasons we expect that in practice the step size will also 
be bounded from below in the obvious manner. This case is not for- 
mally dealt with in the text since the main results may be readily 
inferred from the saturating adaptive quantizer. 

Tor the saturating adaptive quantizer, the following equations 
govern the evolution of {p;(n) = Pr [w(n) = 7]}, 7 S L: 


pilin +1) = Dipepise(n) + ai-upi-i(n) 18 


pin + 1) = a-piiin) L-k+1 81 (13) 


L 
pu(n +1) = ¥ ayp,(n). 





The important super- and sub-martingale properties of the random 
walk, as expressed by the inequalities in eq. (10), apply as well to the 
saturating adaptive quantizer. 


Hi. THE EXISTENCE AND UNIQUENESS OF THE STATIONARY DISTRIBUTION 


We examine in this section questions related to the stochastic 
stability of the adaptive quantizer. We establish theoretically that 
certain acute types of erratic operations such as the unboundedness of 
the evolving random variable, namely, the step size, do not occur. We 
begin by establishing that the process has the basic properties of a 
well-behaved process, namely, irreducibility and recurrence. We 
thereby establish the existence and uniqueness of a finite stationary 
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distribution. We then proceed to the saturating adaptive quantizer, 
the more realistic model of the adapting algorithm, which in addition 
to the above properties, is also aperiodic. Here, the entire state space 
is a single ergodic class. The main result of this section is obtained from 
the construction of a stochastic Liapunov function for the process; and 
the theory of stochastic Liapunov functions is fairly well known.!!-!? 


2.1 Irreducibility of the Markov chain 


The chain is irreducible if and only if every state communicates 
with both the neighboring states. This occurs if and only if there 
exists nonnegative integers m, m’, n, n’ such that 

ml — nk = 1 (14a) 
and 
ml—nvk = —1. (14b) 


It is an elementary fact from number theory that this occurs if and 
only if k and 7 are relatively prime, 1.e., their greatest common divisor 
is unity. In fact, Euclid’s algorithm yields the unknown quantities 
in eq. (14). 


2.2 Recurrence 


Consider the following nonnegative function of the states: 
Vi) = |2| a= 0, +1, ---. (15) 


This function is a stochastic Liapunov function! if the following 
holds: if D(z) is defined as follows, 


ELV{o(n + 1)}|w(n) = 7] — VG) = DO), (16) 


then (7) D(z) is uniformly bounded from above and (72) D(z) S — e€ < 0 
for all but a finite set of states 7. Condition (2) is trivially true for the 
process. Also, for all 7 = k 

; l l 
Di) = ~ +) (%— -45) (+ (bs peq) <0. 


and, for allz < — 1, 


Di) = (k +0) (0: = ro) <(k+1) (0-. . oxi) <0. (18) 


Therefore, condition (72) is verified, and V(z) is a stochastic Liapunov 
function for the process. 
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From IKCushner’s Theorem 7!* we have recurrence* and we can infer 
further, from Theorem 4, that there exists at least one finite invariant 
measure, i.e., stationary distribution. Also, as we have shown earlier 
there does not exist two or more disjoint self-contained subsets of the 
state space; hence, we have from Theorem 5 that there is at most one 
invariant probability measure. Thus, the existence and uniqueness of a 
finite stationary distribution for the step size of the adaptive quantizer 
is established. 


2.3 The saturating adaptive quantizer 


We will circumvent the technical nuisance! posed by periodicity by 
proceeding to the saturating adaptive quantizer. In this case the above 
arguments leading to irreducibility and recurrence are intact. In 
addition, the end state L has period 1 and, since periodicity is a 
class concept (i.e., every state in a particular communicating class 
has the same periodicity), the entire Markov chain is aperiodic. We 
have, then, p(n) — p for any p(0) and p; > 0 for all z. Also, the state 
space is a single ergodic class. Hence, the statistical average of the 
step sizes approach a limit given by the unique, finite, stationary 
distribution. 


iil. SOME PROPERTIES OF THE STATIONARY DISTRIBUTIONS 


In this section we investigate in detail properties of the stationary 
distribution of the step size. In eq. (7) if we set p:(n + 1) = p(n) = pi, 
then the stationary distribution is given by {p;}. Thus, the stationary 
probabilities are the solutions of 


(19) 





> pe =. (20) 


For the saturating adaptive quantizer, we have from eq. (13) that 
the basic recursion in (19) holds for all 7 S (L — k). The remaining 


* A Markov chain is recurrent if and only if every state is recurrent; and state z 
is recurrent if and only if, starting from state 7, the probability of returning to state 
2 after some finite length of time is one. 

Feller’ writes: “The classification into persistent and transient states is funda- 
mental, whereas the classification into periodic and aperiodic states concerns a 
technical detail.” 
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equations are (20) and the following: 
P= a-pit L-k+1sisL-I1 (21a) 


L 
| es zu A;Pj (21b) 


and, of course, p; = 0,2 > L. 


3.1 A useful reduction of the equations for the stationary probabilities 
To provide some insight into the motivation for the step we under- 
take here, consider the recursion, analogous to (19), that would arise 
from a Markov chain with uniform transition probabilities: 
Di = Opin, tapi, atb=1. (22) 


A particular solution of the above recursion is p; = ¢, a constant. Since, 
in probability theory, interest is restricted to solutions with bounded 
sums, one would proceed in the case of (22) by factoring the root at 
unity from the characteristic polynomial: 


bet — Nb + a = 0, 


and thus obtain a new, and reduced, polynomial and an associated 
recursion. This operation is paralleled for the more general recursion 
in (19) by the following: from (19), 


Di — Di-t = OrsnDipn — O¢-1pi-1. 
Hence, for all 7, 


j j 
LX (pi — Pi-1) = Xu (b:44Dit¢~ — 5¢-1ps-1) (22a) 


which reduces to 


(23) 





Remarks: 


(2) Observe that we are justified in carrying out the operation in 
(22a) in the case of solutions of (19) for which >)’. p; is bounded and 
which we have established, in Section II, to be the case for the station- 
ary probabilities. 

(727) The reduction alluded to earlier refers to the fact that the largest 
difference in variable indices in (23) is k + l, while the largest differ- 
ence in (19)isk+1+1. 
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(111) Observe that when k = | = 1, (23) gives the solution in closed 
form: pji1 = (a;/bj41)p; and >¢p; = 1. This is a previously known 
fact; see Feller! and Fine.’ However, neither author gave any indica- 
tion of the possible generalization to the form in (23). 


For the saturating adaptive quantizer, (23) holds for all 7 S (L — k). 
Hence, the range over which (23) is valid is such that every state is 
included in at least one component of the recursion. 


3.2 An identity involving the stationary distribution 


We use eq. (23) to show that the statsonary probability of the nth 
input sample, x(n), falling in the lower slot, Pr, [|x(n)| S A(n)] 
= 1/(k + l). The significance of this identity from the point of view 
of optimal steady-state operation (see Section 1.1) is that by appro- 
priate choice of k and | the above quantity may be matched to the 
corresponding probability for the optimal nonadaptive step size. This, 
of course, has the effect of locating the central step size, eq. (3), close 
to the optimal nonadaptive step size. In the case of independent 
Gaussian inputs, the above quantity is close to 0.68 and a reasonable 
approximation is obtained by making k = 1 and / = 2. 

From (23), 


j+k j 
> bp: = >; Pi- 
j-l+1 jri+l 
Hence, 
00 jak 00 2 
2 ~~ Op = dX > Pi- (24) 
jae i=j—l+1 je—0o i=j—l+1 


The left-hand side equals (k + 1) >>%. b;p;, while the right-hand side 
equals J. Hence, 


(25) 





Consider what the above equality implies in terms of step size 
behavior. The stationary probability of the input falling in the lower 
slot, 


> Pr,[A = Cy and |z| < Cy*] 


t=—=00 


Pr, [ |x(n)| S A(n)] 


DO 0: Pr, [A = Cy] (26) 


~=—% 
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from the independence of {x(n)}. Hence, from (25), 
a e 
k+1 


Immediately on substituting AZ; = y~* and Mz = y' we have an 
identity with a rather appealing and natural interpretation’: 


MyM? = 1 (28) 


where p1 and p2 are respectively the two stationary probabilities of 
the input falling in the lower and upper slots. 
For the saturating adaptive quantizer, it can be shown that 


Pr. [|a(n)| = A(n)] = (27) 


l 
Pm bp, < k+l i (29) 


However, the quantity [(//k + 1) — ¥ bp; ] depends only on (k + 1) 
terms involving the end probabilities pz, ---, pz—z.-1 and it goes to 
zero with these probabilities. Now we will prove in Section 3.3 certain 
results which indicate that these probabilities are relatively small if 
L is large. 


3.3 Geometric bounds on the stationary probabilities 


In this section we prove a fundamental property of the stationary 
distribution of the step size which holds for all values of y. We obtain 
sharp bounds on almost all of the stationary probabilities—the bounds 
apply as well to the saturating adaptive quantizer—which show that. 
the stationary probability of the random walk being in a particular 
state falls off at least geometrically with the distance of that state from 
the 0 state. The actual bounds obtained are substantially stronger and 
they indicate that a localization property on the stationary distribu- 
tion is inherent for the random walk. As discussed in Section 1.1 this 
localization property is important in understanding the basis for the 
satisfactory behavior of the adaptive quantizer. 

We obtain the following point-wise bound: for every 7 > 0 we give 
positive constants 7 > 1 and c such that for all 7 = 2, 


1 \ #7 
pi0(5) 7 (30) 
The quantities r and c depend on 7. The quantity r which we call the 


*D. J. Goodman first conjectured the existence of (28) in the context of the adap- 
tive quantizer. Earlier, N. S. Jayant® made a related conjecture in connection with 
an adaptive delta-modulator. 
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local steepness factor is a monotonic increasing function of 7 for non- 
negative 7. Of course, a corresponding result holds for ¢ < 0 and all 
j 84, 

Let P; denote the (k + 1 — 1)-dimensional column vector* with the 
following components 


P; = [pi Dixt, +++) Ditegi-2)*. (31) 


Then, from (23), we obtain (k +1— 1) X (Kk +1-— 1) transition 
matrices A;, where 
P,., = AP,. (32) 


The leading (& + 1 — 2) components of P;.1 are obtained from P; by 
merely shift operations. The nontrivial information in A; is in the 
last row which is obtained from (23); clearly, A; depends on 7. 

We will show that there exist a constant weight vector 4, every 
element of 2 being positive, and a constant r > 1 depending only on 
A;, such that for all 7 2 7 

wAZ* = rd! (33) 


in the sense that every element of the left vector is not less than the 
corresponding element of the right vector. Since P;;; is a vector with 
nonnegative elements, we have 

TMP 541 Ss AF Py 41 = MP5. (34) 
Hence, 


(35) 





Remarks: Equation (35) is a strong result if 4‘P; is viewed as a norm of 
the vector P; of the Li-type: |x| = >° \;|z:], which is a valid inter- 
pretation since the latter reduces to 4‘x whenever every element of x 
is nonnegative. By standard methods we can obtain upper bounds for 
P; in norms other than the one used in (35). In particular, (30) follows 
trivially. 

It is necessary now to discuss the structure of the matrix A;?. 
Directly from (23) we obtain the first row:t 





(1 — 1) terms k terms 
mii. Ooms —————_—_—__——_r~  ——_— 
_— Gitt Gite —. Gi+l-1 Os41 Ost144 a Ost1ph—1 
a; ? a; } bf a; ) Qa; ? a; P ? a; : 


“The superscript ¢ denotes the transpose. 
T Observe that neither A; nor Aj! is a stochastic matrix (nonnegative elements, 
columns sum to unity). 
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The remaining rows of Aj‘ reflect shift operations: for m = 2, 3, --- 


) 


(AZ))mn = 0 if n¥ (m—1) 
=] if n=(m-—1). 


Before proceeding to prove (338) we need the following lemma. 


This lemma concerns the matrix Aj! which is obtained from Ap?! by 
merely replacing the first (1 — 1) elements of the first row by —1. 


Lemma 1: For every 1 = 0 


(i) Aj} has a unique positive real eigenvalue r, say. Furthermore, 
n> 1: 

(21) Every element of the corresponding left eigenvector 2% is of 
the same sign and nonzero; hence, 2 may be taken to be a 
positive vector. 

(222) r, which depends on 72, is monotonic, strictly increasing with 7. 


We give the proof of Lemma 1 in Appendix A. 
We need one further observation to prove (33) with the help of the 
lemma. For 7 = 1, 


AZ? = a4(AZ! — Az!) + ata? 
= x4(Az! — Ap) + rat. 
The bound in (33) follows if 1“(A7! — Aj!) = 0. Since 4 is a positive 
vector it is sufficient to show that the elements of the matrix 


(Az! ~ Az") are nonnegative. The only nonzero elements of the 
j 


matrix (Aj? — Aj?) are in the first row. That every term of the first 
row is nonnegative is implied by the following: for s 2 1 





jee SG (36) 
Qj 
and 
bite _ Dita x g, (37) 
Qj a; 


This concludes the proof of (83) and, hence, of (35). 


Remarks: 


(2) The reader may now appreciate the reason for replacing some of 


the elements of Aj! by —1 to form Aj’: a,,;/a; although bounded by 
1 can come arbitrarily close to 1. 
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The reader is also due an explanation for our having worked with 
Aj’ after defining the natural transformation Aj, especially since 
(34) may be put in the form A‘{I — rA; JP; 2 0. The reason is that 
r and i, depending only on 7, do not exist such that for 7 = 7, 
X41 — rA;|20, although, as we have shown, 4 and r do exist such 
that a‘LI — rA;]A;7' 2 0. In working this step the assumption of 
P;,; 2 0, rather than P; 2 0, appears to be critical. 

(22) The interesting quantity r = r(z) may reasonably be called the 
local steepness factor, since for 2 2 O it is a local measure of the rapid- 
ness with which the stationary distribution falls off. From statement 
(272) of the lemma we have the fact that the distribution tends to get 
steeper with increasing distances from the natural center of the dis- 
tribution, the 0 state. 

(177) The theoretical interest in the inequality in (35) results from 
the fact that we cannot expect to obtain a significantly better value 
than r for the geometric factor in geometrical bounds on p; for all 
j = t. The reason for this is that by making 6,41 very close to b; over a 
fairly large set of 7’s, it is possible to make the solution of (23) close to 
the stationary probabilities of a random walk with uniform transition 
probabilities, which in turn may be obtained in terms of r as the unique 
positive real root of the characteristic polynomial C(u) given in 
eq. (56), Appendix A. 

(iv) From symmetry we expect results similar to (35) to hold for 
7 < 0. Perhaps the simplest way to show this is by means of the follow- 
ing transformations which have the effect of making the direction of 
decreasing 7 the forward direction. Let 


, ; 
P+ = Dy, 0-4~-= 4; a~4~= D;. 


The basic recursion (23), stated in terms of the new variables, is 
j+l , 4 j / ? 
Dw bp = DL C1 — 0i)pi. 
iI j-k+1 


Now {b;} is a monotonic, increasing sequence with 7 and i > 0 = Ib; 
> ka;. (Observe the interchange of | and k, ie., Il’ = k and k’ = 1.) 
This transformation makes the transfer of results holding for 7 > 0 
to z < 0 fairly straightforward. 

(v) In considering the application of (35) to the saturating adaptive 
quantizer we note that the basic recursion (23) holds over the entire 
range of states, i.e., (23) holds for all 7 S L — k. Hence, (35) holds for 
L-—-(U+k)+22j2720. This observation is the basis for a 
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statement made earlier in Section 3.2, namely, we expect the tail 
probabilities of the stationary distribution of the step size for the 
saturating adaptive quantizer to be small. 

From (35) we obtain a rather simple point-wise bound on the 
stationary probabilities. Let \,, denote the largest element of the vector 
x. Clearly, * 

APES Anl'P:, 


and, hence, from (35), for all 7 2 7 2 0 


: ge% 
Km Ss a P; < ( WP; Ss (=) Am(1‘P;), 


1.€., 


Hence, 


J I 
< ] (38) 





where r = r(t). 


3.4 Lower bounds on the steepness factors, r(i) 


We have associated with every state 7 a local steepness factor r(z). 
Here we go back to the definition of r(z) as being the unique positive 
root of the polynomial C(u), eq. (56), to obtain the following bound 
which has the advantage of being explicit. 





. 1/(k+1—1) 
a | A210, FeO: (39) 


la; 


Observe that p(z) > 1 for all 7 >0 and itself forms a monotonic 
increasing sequence with 7. To prove (39) it is enough to show that 
C((kb;/la;) |] S$ 0. The proof is straightforward but tedious and we 
omit it. 


3.5 The effect of y on the stationary distribution 


We show here that the mass of the stationary distribution of the 
step size can be localized about the central step size to an arbitrary 
extent by making y sufficiently close to unity. To do this, we first put 


*The column vector with every element equal to unity is denoted by 1. 
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together from the results of the preceding sections a rather explicit 
bound on the stationary probability of the step size exceeding a par- 
ticular value for a given y, i.e., Pr; [A = Cy*]. This bound is in a form 
which allows direct comparison with the corresponding probability 


arising from the choice of y’ = Vy. By successively taking y to be 
the square root of the preceding value, the bound on the probability 
can be made as small as desired. As before, we shall restrict our at- 
tention to step sizes which exceed the central step size, i.e., 2 > O since 
a parallel argument holds for z < 0. 

For 7 > 0 and r = r(t), we have from (35) that 














(Dr) dD ps ¥ WP; S A'P; ¥ (;) ‘= - (40) 
j=it+k+l—2 j= j=0\ 7 1 
Now, as in (39), 
an [ fbi) MED 
re p(t) = ( a) 
and 
xP; 
DA, S max [pi +++, Dizeti2]. 
Since 
Pr.[A 2 CyH2] = Yop, 
| j=itk+—-2 
we have, from (40), 
Pre[a 2 Cy##t2] = 2 max [ps +++, pesasral | (41) 
Finally, from (88), fora 2k +1 —1, 
Se ene ees ea Ges (42) 
| a) » Pitk+t— = p(1) . 





Equations (41) and (42) give the bound for the mass of the distri- 
bution to the right of a particular state, which we shall now compare 
with a similar bound that holds for y’ = Vy. The prime superscript 
will be used on symbols to denote the functional dependence of the 
associated quantities on y’. In establishing the reference (central) step 
size [see eq. (3) ], minor differences exist depending on whether 


@) Pr[|x(n)| Sy] < = Pr[|z(n)| = y*"?] 


ris 
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or 


Gi) PrCle(n)| S ¥] <5 S PrLle@)| S74) 
We consider only (22), in which case: w’(n) = 245 w(n) = 7 and 
bs, = 0; for all 7 = O. 

Repeating the arguments leading to (41) and (42) we have 


; ? 2 D , 
Pr, [A = Cy y2it e+1—2)7 < voy I max [ pei, oes, Poi+k+1—2 | (43) 


and 
; ‘ 2 1 2i—k—l 44) 
max [ poi, seer Poi+e+i—2 | = 2'(2) | : ( 
Since p’(2z) = p(t), we have 


i+ (k+1—2 p(t) > t—k-l+1 sake i-l 
Prefaz ovens | | Lao] 


Comparison with (41) and (42) completes the demonstration. 


IV. TRANSIENT RESPONSE 


The preceding section discusses various aspects of the stationary 
distribution of the step size which effectively describes the steady-state 
behavior of the device. However, as stated before in Section I, the 
steady-state response is only of partial interest since the adaptability 
of the device is tied to quickness of response in the following situations: 


(2) Start up—we are forced to consider situations in which the 
initial step size is fairly arbitrary. 

(72) Changes in the scaling of the input distribution—the scenario 
here is that the device has adapted to a particular intensity level 
(scaling) of the input distribution when a jump occurs to a new 
intensity level. 


In common with both situations, we have an initial step size and a 
waiting time for the step size to adapt to the desired step size. Recall 
that with k and / appropriately chosen, the desirable step size is the 
central step size, which corresponds to the 0 state in the random walk, 
eq. (5). This aspect of the behavior of the device is also related to the 
rate at which the evolving step size distribution approaches the station- 
ary distribution. 

The main contribution of this section is the development of formulas 
for the efficient computation of the mean time required for the step 
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size to first reach the central step size for various values of the initial 
step size. The designer can use the information generated by the 
methods given here in the following manner. Assuming that the de- 
signer has some understanding of the rate of variation of the intensity 
level of the input distribution, he is in a position to determine the 
smallest value of y for which the adaptation algorithm adequately 
tracks the input process. The parameter y has to be made sufficiently 
large for the mean waiting time (time, of course, is used synonymously 
with number of transitions) for adaptation to be small compared to 
the changes in the location of the desired step size arising from changes 
in the intensity level. 


4.1 The mean time for first passage to the origin 


We will consider the random walk, eqs. (5) and (12), for the satur- 
ating adaptive quantizer since in the limit, as L becomes large, the 
functionals obtained for this model yield corresponding quantities for 
the adaptive quantizer. Also, we shall consider only the case of the 
initial state w(0) > O since the results obtained can be transferred to 
the case of negative initial states in a fairly obvious manner (see 
Remark (iv) of Section 3.3). 

Let the initial state w(0) = 7 > O and let M; denote the mean time 
required for the first occurrence of the event w(n) S 0. We observe 
that for all values of L, not necessarily finite, the time to first passage 
is finite with probability 1 as a consequence of the properties of recur- 
rence and irreducibility established earlier in Section II. If the first 
transition results in a decrease of the step size, the process continues 
as if the initial state has been (2 — k). The conditional expectation of 
the first passage time, therefore, is M/;-, + 1. From this argument we 
deduce that the mean first passage time satisfies the recursion* 





The relation in (46) may be used to generate the entire sequence {M,} 
provided the initial conditions are known. Now, by the same argument 
that led to (46), we have that (46) holds for1 S$ 7 S k with 


Mi+.= Mo,» = --: = Mo = 0. 


“There is some similarity between (46) and the equations arising in gambler’s 
ruin problems! and sequential analysis,!* in which generally k = | = 1 and the transi- 
tion probabilities are not variable. 
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The remaining / boundary conditions, namely, 
Me Mig tos, MG 


are hard to obtain and it is necessary to look more deeply into the 
dynamics of the process to obtain these quantities. 

For every sampling instant we define the L-dimensional vector z(n) 
with components z;(n), 1 S 7 S L, where 


z;(n) = Pr[w(n) = 7 and w(s)=1 foralls < nj. (47) 
These vectors, z(n), evolve with time according to 
zin+1)=Dz(n), n2 0. (48) 


These equations are given in Appendix B. Here we reproduce the 
structure of the L X L matrix D: 


k 


SSS 
Oy ths 0. Bea 
i+ |: Dy42 


0 
a1 ; 
D = as * 
: L 
0 
0 
Qp-l QAn-l41 °'* Gy 


Putting together various properties of the matrix D and the random 
walk, we obtain, in Appendix B, the following result: for 1 = 1 


(49) 


where [I — D]x® = e; 





and the elements of the vector e; are zero everywhere except at the 
ith location where it is unity. In Appendix B it is shown that [I — D] 
is nonsingular. We observe parenthetically the virtue of the recursion 
given in (46) in that it allows us to generate rather easily all the 1/,’s 
once the / inversions necessary to evaluate M,, ---, Mi are carried out. 

The matrix inversion in (49) may be viewed as a mixed boundary 
value problem with the first 1 and the final & equations providing the 
boundary conditions. The bulk of the elements of the vector x“ 
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satisfy a recursion that was encountered previously in Section IIT: 
xy = dipetf te + a;—ejle. (50) 


Furthermore, we show in Appendix B that the elements 2% are all 
nonnegative. Hence, we are in a position to usefully apply, even for 
infinite L, the techniques and results of Section ITI. 

First, we carry out the reduction of the equations as stated in 
Section 3.1 where the motivation for this step is discussed. We obtain 


r+k r 
25 Oia =... 25. “C= b)ay bar S Ga bh), (51) 
j=rt+l1 j=er—l4+1 


The superscripts on the x’s have not been used since (51) holds for 
alx®, I SiS. 

One benefit of the above form is that it involves one less variable 
than the original recursion (50). In the important case of k = 1 and 
l= 1, this reduction is sufficient to transform the original mixed 
boundary value problem (49) to an initial value problem, i.e., the solu- 
tion to the matrix inversion problem (49) satisfies a recursion with 
specified initial conditions. Exact computation in this case becomes 
quite trivial. The details of this solution are given in Appendix C. 
Apart from its independent interest, this result is of particular interest 
in the adaptive quantizer when the distribution of the input sequence 
is Gaussian. As discussed previously, it is desirable to have in this 
case 1/(k + 1) = 0.68, and k = 1 andl = 2 will suffice. 

Another property of the solutions x“ of (49) which holds for all L 
is that with increasing 7, x} decreases at least geometrically. This con- 
clusion may be drawn from the bounds obtained in Section 3.3, eqs. 
(85) and (88). From the point of view of numerical inversion of 
[I — D] for large L, this is a critical property in that it is a necessary 
condition for most numerical techniques. The reader is referred to 
Richtmyer and Morton!’ for one such technique that we have used 
successfully and found to be efficient in that it effectively exploits the 
band structure of the matrix [I — D ]. 

Finally, we remark that while we have dealt exclusively with first 
passage across the 0 state it is clear that generalizations to first cross- 
ings across states other than the 0 state is straightforward. 


4.2 Bound on the mean first passage time 


Two formulas, eqs. (46) and (49), have been given for computing 
the mean time required for the step size to adapt from an arbitrary 
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initial value to the desired, and also central, step size. However, by 
examining these formulas it is not easy to gain insights into the rate 
at which this adaptation time grows with the distance separating the 
two states and its dependence on y. Here, by probabilistic reasoning, we 
obtain an explicit upper bound on this time and this bound does pro- 
vide some insight. As we have done before, we consider here only the 
case of positive initial states, 1.e., w(0) > 0. Let M;;,0 S 7 < j, denote 
the mean first passage time under the following conditions: the initial 
state w(0) = 7 and first crossing occurs after 7 transitions if w(r) S 7 
and w(n) > 7 for all n < 7; then M,; = E(r). In this notation the 
quantity 7; defined in Section 4.1 is equivalent to Mo,. 
In Section 1.2, eq. (10), it is given that 


BGA Aaa e. 2 rt (52) 


Denote the quantity on the right by —S; and observe that for 7 > 0, 
Siz1 > S: > 0; hence, the supermartingale property. [For the saturating 
adaptive quantizer, the supermartingale property holds even more 
strongly, i.e., for 7 > 0, (52) holds with the equality replaced by S. ] 
In fact, the supermartingale property holds for the transformed 
process: w’(n) = w(n) + nSi41. Le., 


E[w'(n + 1)|w’(n)] S w’(n) (53) 


for all w’(n) = (@ + 1) + nSiii. For the crossing problem, (53) holds 
for all (n + 1) S 7, the crossing time. We can now apply a theorem 
due to Doob?8 on optional stopping on supermartingales. In this case, 
the theorem states that 


E[w!(r)] S$ Elw"(0)]. (54) 
Since 
G@+1—b) + Sirk(r) S Elw'(s)] S$ Elw'(0)] = j, 


we obtain 


(55) 





Mj; = E(r) < gal Hie) y 


We gain some insight on the role of y in determining the transient 
response of the device by observing the dependence of the above bound 
on y. Suppose we are interested in Mo,, the waiting time for the initial 
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step size A(0) = Cy’ to reach the central step size C’. Consider the 
effects of making y’ = vy on this waiting time (the multiplier coeffici- 
ents of the device are therefore Vy-* and vVy!). We let the prime 
superscript on symbols indicate a functional dependence on y’. In 
establishing the new central step size [see eq. (3) ], minor differences 
exist depending on whether 


A 


@) Priam)! Sy) < po S PrCle@)| S74] 


or 


(1) Pr[|x(n)| 


IA 


l 


ras Priam) 


y*]. 


IA 
IIA 


yril< 


We consider only (72), in which case the central step sizes are identical: 
w'(n) = 22 w(n) = 7 and by; = b; for all ¢ = 0. The waiting time 


y = 1.11/2 
M, =0.95 
My a 1.1 


MEAN TIME FOR FIRST ARRIVAL OF STEP SIZE TO CENTER 





INITIAL STEP SIZE 


Fig. 3—Transient response of the adaptive quantizer. 
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y= 1.1 1/2 
M, = 0.95 
Mo = 1.1 


y= Al 
M, = 0.909 


MEAN TIME FOR FIRST ARRIVAL OF STEP SIZE TO CENTER 


0.2 0.3 
LOG yg (INITIAL STEP SIZE) 





Fig. 4—Transient response of the adaptive quantizer. 


for the step size to adapt from identical initial step size Cy to final 
step size C is Mo.2;. From (55), 


, Ptaaet 
Moe; gil2d — (k — 1)]. 
1 


changed has the effect of making the bound on the waiting time at 
least twice as large for 7 > k. This is a conclusion which is plausible 
in the light of the linear form of the bound (55) since the effect of 


making y’ = vy is to introduce twice as many transitions between the 
initial and final step sizes. 


Now, So S Si S Si; hence, making y’ = Vy and keeping k and 1 un- 


4.3 Computational results 


We present here a sampling of our computational results. It is 
assumed that for every n, x(n) is normally distributed with unit 
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variance. The optimal step size A in this case has the property that 
Pr {|x(n)| S A} = 0.68. To center the stationary distribution of the 
step size close to the optimal step size, we choose k = 1 and l = 2. 
Figure 3 plots the mean time for first passage to the optimal step 
size vs. initial step size, and the initial step sizes chosen for this figure 
exceed the optimal step size. Various values of y(M1 = y~*, M2 = y') 
were used. Figure 4 provides the same information except that the 
horizontal axis corresponds to logiy A(0), rather than A(0) as in Fig. 3. 
The mean first passage times 7, and M, were obtained by the method 
outlined in Appendix C, and M;, 7 = 3 were generated by using the 
recursion in (46). To give some idea of the rate of convergence for 


a, eqs. (70) and (71), we tabulate some values of 2x} for the case of 


MEAN TIME FOR FIRST ARRIVAL OF STEP SIZE TO CENTER 





0.5 
INITIAL STEP SIZE 


Fig. 5—Transient response of the adaptive quantizer. 
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y = 1.1: 


ge [ck .2 3 4 5 6 7 8 9 10 
2f: | 1.4 0.53 0.66 0.381 0.20 0.08 0.03 0.92 K 107? 0.24 & 107? 0.41 X 1073 





ji 11 12 13 14 15 16 
xf: | 0.59 X 10-4 0.53 & 10-° 0.35 K 10-® 0.13 K 1077 0.30 X 10° 0.31 K 10-4 
Figure 5 is similar to Fig. 3 except that here the initial step sizes are 
less than the optimal step size. Figure 6 plots the same information 
with logio A(O), rather than A(O), on the horizontal axis. The mean 
first passage time M, was obtained by solving (49) by the method 
given in Ref. 17 and all other first passage times were generated by the 
recursion in (46). 


MEAN TIME FOR FIRST ARRIVAL OF STEP SIZE TO CENTER 





-0.6 -0.4 
LOG qo (INITIAL STEP SIZE) 


Fig. 6—Transient response of the adaptive quantizer. 
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APPENDIX A 
Proof of Lemma 1 
Proof: 
(i) Az} being in the form of a companion matrix, the coefficients of 


the characteristic polynomial of the matrix are the elements of the 
first row: 





C(u) & (—1)# det [AP = uD] 
= prtl-l it... + yk — [ aye} + agp? + ++) + ax |, (56) 
where 
ae Dist Br Dietyy ee) aR = Diptek—t | (57) 
a; a; Qi 


By Descartes’s rule the polynomial C(y) has at most one positive real 
root. Since C(0) = — a, < 0 and C(un) > © as p— ~, there exists 
exactly one positive root. Let r denote this root. 

Now C(1) < 0 if la; < (b541 + Des t41 -}- ae + brte—1)- The latter 
condition holds for all z = 0. Hence, r > 1. 

(it) The left eigenvector 4 corresponding to the eigenvalue 7 satisfies, 
by definition, afAZ? = rv. Examining the component equations we 
find that 


VE eee ee ye) Pere. (58) 
Also, 
Nes | 
hipk-t = mabe int Ae pee Ae ee ae ve) AS Ske 159) 


Finally, r\1:.-1 = adi. Since the a’s are positive quantities, the state- 
ment is clearly true. 
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(22) The statement can be verified by inspecting the characteristic 
polynomial C(u) and using the fact that the coefficients aj, ---, a, 
each increase with 7. 


APPENDIX B 
Derivation of equations (48) and (49) 


The derivation of the equations governing the evolution of the 
vectors z(n) defined in eq. (47) proceeds as follows. For convenience, 
let X(n) denote the event 1 S w(r) S L for all 7, 0 S$ 7 S n. Hence, 
by definition, 


zn) = Prf[w(n) = 7 and X,] I1SJISL. 
Since 


z;(n) 


Prfw(n) = 7 and Xy_-1] 


I 


© Pr Co(n) = jlo(n — 1) = 6, Xeale(n — 1) 


bias Zeag(n — 1), 1SEJS1, 


a;-1 2j-1(n — 1) + Opzye(n—1), C+DS78S (Lk) 
= <4;-1 z;-(n - 1); (L ee 1) Ss j S (L ie 1), 
ZL 
>» aaein—-—1) jg=L. 
<r 


The above equations define the matrix D which relates z(n) to z(n — 1) 
as in eq. (48). 

For the derivation of eq. (49) we proceed as follows. For 2 = 1, 2, 
-++, DL, let 


F(n + 1) & Pr [first passage occurs at (n + 1)|w(0) = 7] 
Prfw(n+1) S50, X,|w(0) = 7] 


I 


k 
= >) b;2,(n) with z(0) = e. (60) 
j=1 


The vector e; has every element equal to zero except for the 7th element 


which is unity. To express eq. (60) in vector form we let b = [bybo-+ by 
0---0]!. Then, from (60), 


Fj(n+ 1) = b'z(n) with z(0) = e:. 
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By definition, we have that the mean first passage time conditional on 
the initial state being 2, 


M,; = 2 (n + 1)F,(n + 1) 


bé 2 (n + 1)z(n) 


= bt >> nz(n) + dt 2, 2(n). (61) 


n=0 


Now the second term in the above expression is unity since the proba- 
bility that passage occurs at finite time 1s unity. Now consider 


fI—D] a nz(n) x nz(n) — 2 nan + 1) 


I 


a z(n) — z(0). (62) 


Hence, denoting by 1 the column vector with every element equal to 
unity, we have from (62) that 


1‘{I —D] a nz(n) = I! 2 z(n) — 1 (63) 
= bt a nz(n), (64) 


since 1'z(0) = 1 and bt = 14[I — D]. It only remains to consider 


>», zZ(n) = pa D‘| z(0). 
n=0 i=0 

The above series converges since every eigenvalue of the matrix D 
lies strictly within the unit circle in the complex plane. The proof of 
this follows from an old matrix theorem!® which states that if the 
diagonal elements of the columns weakly dominate the sum of the 
absolute values of the off-diagonal elements with strong dominance 
holding for at least one column and the matrix is irreducible, then the 
determinant is nonzero. Applying this theorem to [D — AI], |A| = 1, 
we note that the irreducibility of the original Markov chain implies 
irreducibility of the matrix [D — AI] and that the weak column 
dominance property holds everywhere while the strong column 
dominance property holds for the first k columns. Hence, 


2 2(n) = x D*| z(0) = [I — D }“'z(0). (65) 
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Putting together the above results we have (49), namely, 
M;= > «x where [I — D]x® = e,. 
jz1 


Observe that x‘ = Yz(n) and, from the definition of z(n), it follows 
that every element of x™ is nonnegative. 


APPENDIX C 
Mean first passage times for the case k=1, |=1 


We have as our starting point eq. (49), namely, 


M;= 2s a, (66) 
j 


where [I — D ]x® =e, (67) 


and we are interested only in 1 S71. 

The transformation that was made in Section 3.1 is equivalent to 
the following: add to each row, 7, of [I — D Jallrowsr + 1,r + 2,---; 
and do the same to the vector e;. This operation makes the matrix 
[I — D ] lower triangular, the reason being that with the exception of 
the first column, the elements of all other columns of [I — D] sum to 
zero. The resulting equations are as follows: the first component 
equation yields 

by? = 1, (68) 


and the next (J — 1) equations: 2 Sr S 1, 
1 if ri 
r—1 
a ee (69) 
j=l 
0 if r>i. 
Finally, 


r—l 
gh = 2 Ds age fore ST. (70) 


b, j=r-—l 





The boundary conditions to the basic recursion in (70) are in (68) and 
(69) which are, of course, solvable: 







—_— 
IA 
3 


* 
<1 @ = 1/ [J }; 
j=l 


: (71) 
sto oa? =(@P—1)/ IL 53. 
j=ttl 
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Spectra of Digital Phase Modulation 
by Matrix Methods 


By V. K. PRABHU and H. E. ROWE 
(Manuscript received October 31, 1973) 


We derive the spectral density of a sinusoidal carrier phase modulated 
by a random baseband pulse train in which the signaling pulse duration 
us finite and the signaling pulses may have different shapes. The spectral 
density 1s expressed as a compact Hermitian form in which the Hermitian 
matriz 1s a function of only the symbol probability distribution, and the 
associated column vector 1s a function of only the signal pulse shapes. If 
the baseband pulse duration is longer than one signaling interval, we 
assume that the symbols transmitted during different time slots are 
statistically independent. The applicability of the method to compute the 
spectral density is illustrated by examples for binary, quaternary, 
octonary, and 16-ary PSK systems with different pulse overlap. Similar 
methods yield the spectral density of the output of a nonlinear device 
whose input 1s a random baseband pulse train with overlapping pulses. 


l. INTRODUCTION 


In recent years, digital phase-modulation techniques have been 
playing an increasingly important role in the transmission of infor- 
mation in radio and waveguide systems. Various methods have been 
developed recently for computing spectra of a sinusoidal carrier phase 
modulated by a random baseband pulse train.!~§ 

In this paper, we derive the spectral density of a carrier phase modu- 
lated by a random baseband pulse stream in which the signaling pulse 
duration is finite and the signaling pulses may have different shapes. 
The spectral density is expressed as a compact Hermitian form in 
which the Hermitian matrix is a function of only the symbol proba- 
bility distribution and the associated column vector is a function of 
only the signal pulse shapes. If the baseband pulse duration is longer 
than one signaling interval and the pulses from different time slots 
overlap, we assume that the symbols transmitted during different 
time slots are statistically independent. 
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The present method also yields the spectral density of the output 
of a nonlinear device whose input is a similar baseband pulse train. 

The work reported here generalizes and simplifies prior results. The 
form of the present results provides an appropriate division between 
analysis and machine computation that enhances physical under- 
standing and simplifies numerical computations. We compute the 
spectra of binary, quaternary, octonary, and 16-ary PSK systems, 
with overlapping baseband modulation pulses of several shapes. 


If. M-ARY PHASE-MODULATED SIGNALS 
We seek the spectrum of the digital phase-modulated wave 


a(t) = cos[2rft+ ¢(t)], f. > 0, (1) 


where 


60 = YX gult— kT), se = 1,2, +5 M. (2) 


The discrete random process s;, is assumed strictly stationary; as 
noted in (2), it takes on only integer values from 1 to MM. The carrier 
frequency is f.. The signaling alphabet consists of time functions, 
91, 92, °**, gm, that may have different shapes; one of these is trans- 
mitted for each signaling interval 7, to generate the digital baseband 
phase modulation ¢(¢). The different signaling waveforms in (2) may 
overlap, and may be statistically dependent throughout the present 
section and Appendix A. 
For convenience, we define?® 


v(t) =e79(t); (3) 
then 
a(t) = Re {e%?*fe! y(t) }. (4) 


Appendix A shows that the spectrum of x(¢) is 


P.(f) = aP(f — fe) + aP (— f= Jo); 
Qf. # ms MS Ay 2y 35", (5) 


where? 
Pf) = [ ba(neverdr, (6) 
a 1 74 
&,(7) = 6,7) = lim =| ; &,(t, r)dt, (7)" 
A~ 2A —A 
©, 7) = (0 + zr) v*(t)) = (etleGrn-eel). (8) 
t The symbol denotes average on ¢ throughout. 
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The first term of (5) is the spectrum of the complex baseband wave 
v(t) shifted to the carrier frequency + f.; the second term is the spec- 
trum of v*(¢) shifted to — f,. The spectral relationship of (5) is strictly 
true as long as the condition on f, is satisfied, whether or not the two 
spectral terms overlap; that is, this result applies strictly to both 
narrow-band [f, >> bandwidth of P,(f)] and wideband modulated 
waves. 

The condition of (5) requires that twice the carrier frequency is not 
an integral multiple of the signaling rate 1/7. P,(f) has in general 
both continuous and line spectra, the lines occurring at frequencies 
n/T for integer n. The condition in (5) guarantees that the line com- 
ponents of the two spectral terms never coincide. In the exceptional 
case, where the two sets of lines do coincide, it is not surprising that 
it no longer suffices merely to add powers of the two terms, as in (5); 
however, if f, is high enough the modulated wave is narrow band, the 
two spectral terms of (5) do not overlap significantly, and (5) pro- 
vides a good approximation even if the condition is violated. 

Note that we have not found it necessary to randomize the phase 
of the unmodulated carrier or the position of the time slots of the 
digital modulation, as is done in various other studies. Thus, rather 
than (1) we might have considered 


x(t) = cos [2rf + o(t — to) + do], (9)" 


with ¢(¢) still given by (2). The spectral relation of (5) will again be 
strictly true when.2f,7' = integer if ¢o and f are independent, and 
either is uniformly distributed over a suitable interval (Appendix A). 
However, we retain the original formulation of (1) and (2) throughout 
—corresponding to ¢o = 0, tp = 0 in (9)—to determine the effect of 
deterministic carrier and modulation phase on the statistics of the 
modulated wave. 

Consequently, we study only P,(f) throughout the remainder of 
this paper. This suffices for all cases except a low carrier frequency f, 
that is a precise multiple of half the baud rate 1/27; the additional 
calculations necessary for this exceptional case are suggested in Ap- 
pendix A, but not carried out in detail. 


Ill. NOTATION AND STATISTICAL MODEL 


The introduction of vector-matrix notation greatly simplifies the 
ensuing analysis. 


' Hither of the two parameters, ¢o or to, shifts the relative phase between carrier 
and modulation; thus, only one of them is really necessary. 
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First, we rewrite the phase modulation of (2) as 


o(t) = z fafgitt — kT) + af get — kT) 
+--+ + af%u(t — kT)]. (10) 


For a given k (1.e., for a given time slot), one of the a,’s is unity and 
the rest are zero; 
a{s® a 1, 


(11) 


where s; is the strictly stationary, discrete random process of (2), 
taking on the integer values 1, 2, ---, M. 
Now we define for convenience the row vectors 


an, = Lar? af? +++ a’? J, 
et) = Ln ge) +--+ gard) J, 


whose components are respectively the coefficients and pulse shapes of 
(10). The transposes? of these row vectors are the column vectors 


(12) 


ay” gi(t) 

; qg® ; t 
aj=a'={|% |, e@]=e0,' = |% |. (13) 

aft? out) 


Define the unit basis row vectors 
e,=(1 00 ---0, 


,=[0 10 --- 0], (14) 


eu, = L0 O--- O 1], 


with corresponding transpose column vectors e; |, e2 |, ---, ea]. Then 
ax, takes on only the values @€1, €2, --+, €ar, by (11): 


an = lay (15) 


Now we rewrite (10) in vector notation as 


6 =X angle — k7)], (16) 
where - signifies ordinary matrix multiplication throughout. The term 


t Boldface quantities denote matrices throughout. Row and column vectors are 
distinguished by the additional notation -4 and ], respectively. 
* The transpose of a matrix is indicated by the symbol ’. 
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ax,is a vector-valued, discrete’ random process, strictly stationary by 
the assumed strict stationarity of s, of (2). We define the first- and 
second-order probabilities of ax,or ax_] as 


w; = Pr{a, = e;} = Pr{s, = 7} = 0. (17)! 
W(t, J) = Pr{an = Ci, Angn = €;} = Pr{sy = 1, Stun = J} 2 O. (18)? 


w; is the probability that the 7th signaling waveform g; is transmitted 
in any time slot, W,(2, 7) the joint probability that the signaling 
waveforms g; and g; are transmitted in two time slots separated by n 
signaling intervals. w;and W,(2, 7) are independent of k by the assump- 
tion of stationarity. Then since the marginal probabilities are obtained 
by summing over the joint probability function, 


M M 
X W(t, J) = Wj; 2d W(t, )) = Wi. (19) 
i= j= 
Normalization of the total probability to 1 requires 
Mf Ma M 3 
ey wes, 2d dX WG, 7) = 1. (20) 
i= i=1 j= 


Now we introduce vector-matrix notation for the probabilities. Let 
Ww, = [ wy Wg. t8* wa | (21) 


be the probability row vector whose elements give the probabilities of 
the different signaling waveforms, with transpose column vector 


w]=w,’. (22) 
Let 
W.C,1) WaGl,2) +--+ Wal, df) 
W (2,1) W,(2,2) ++: W,(2,M) 
W, = : ; bee (23) 


be the matrix whose elements specify the joint probabilities of all 
pairs of signaling waveforms separated by n time slots. Further, define 
Heil =f1 1---1] (24) 


as a vector with all 17 elements unity. Then (19) and (20) may be 
written as 
LW,=w W,-1l) = w] (25) 


T az,is defined only for integer values of its independent variable k. 
tIt is understood that a, e are either row or column vectors throughout (17) and 
(18), and in succeeding equations. 
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and 
I-wl=wij=1, 1:W,:1J=1. (26) 


The probability matrices have the following useful properties. For 
n = 0, the joint probability matrix of (23) and (18) is diagonal, with 
diagonal elements the first-order symbol probabilities; 

Wy 0 
Wo= fo = Wa, (27) 
0 WM 


where we define the diagonal matrix wa for later convenience. Then 


wal]=w], L:wa = mM, (28) 
By symmetry, 
W.(2, 7) = W-2C9, 2), (29) 
or in matrix notation 
Wi Ws. (30) 


We assume that signaling waveforms in widely separated time slots 
are statistically independent: 


lim W,(4, 7) = w; U3, (31) 
or in matrix form 
lim W, = w]-w, (32) 
since 
(w]-w)’ = w]-w, (33) 
(80) and (32) yield 
W,= W. = W_,. (34) 


Using the above, the mean and covariance of the strictly stationary, 
vector-valued, discrete random process a; are 


(ax ]) = wl, (ax) = My, (35) 
@,(n) = (aryn “an) = W,,. (36) 


The assumption of independence in (31) and (82) implies uncorrelation 
in (86) as n—o. Conversely, uncorrelation implies independence, 
because the random vectors az are unit basis vectors (15),1° rendering 
the covariance and probability matrices identical. Thus, rather than 
assuming independence in (31) and (32), we might equally well 
assume that the modulation vectors a, become uncorrelated for 
widely separated time slots. 
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The spectral density of a; is the discrete transform of (36); 


Piaf) = dX ePrsow,,. (37) 
Equations (36) and (37) are the matrix extensions for discrete vector 
random processes of similar scalar relations for discrete scalar random 
processes;!! the symbol ~ indicates that a discrete random process is 
under consideration. We separate (37) into line and continuous spectral 
components: 


P,(/) = P.,(/) a P..,(/f). (38) 

Then from (32) and (87): 
Pu(f) = ww, X af»), (39) 
P.(f) = DPW, — wm). (40) 


The assumption of (31) and (82) that signaling waveforms in widely 
separated time slots are independent, or the equivalent assumption 
that az,, and a, become uncorrelated for large n, eliminates line 
components except de if we confine our attention to the fundamental 
frequency interval |f| < 4, and so retain only the n = 0 term in 
(39).* Consequently, the modulation has no periodic patterns. 


IV. PSK AS A BASEBAND PULSE TRAIN 


Our problem has been reduced to determining the spectral density 


of the complex wave v(t): 


6) = 2X aye — kr], (42) 


where a and g are M/-dimensional vectors. We show that if the signal- 
ing pulses are strictly time-limited to an interval KT’, v(t) may be 
written 


ot) = DX bust(t — kT). (43) 


' While this relation is defined for the entire range — » < f <~, P, (f) is periodic 
in f with unit period, and only the fundamental period | f |< 3is normally of interest. 
Thus, the inverse transform 1s 


$,(n) = Wa = [ B, (fetrmay. 


+ We may thus write a, = w + acz with (40) giving the spectral density of a. and 
(39) the spectral density of w. 
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Vectors b and r are M*-dimensional vectors expressed in terms of a 
and g, respectively. The different terms in (43) are strictly nonover- 
lapping: 

r(i)]=0], ¢50, t>T. (44)t 


Each b; is a unit basis vector, i.e., only one of its M@* components is 
unity, all others being zero. The spectral density of (43) is determined in 
Appendix B. 

Figure 1 shows portions of ¢(t) of (42) for four different maximum 
signal pulse durations; the terms k = — 1, 0, 1, 2 of (42) are shown, 
and for convenience a; has been taken the same for each of these time 
slots. The pulses are positioned along the time slots such that the 
limits of each signal pulse lie on the boundary between adjacent time 
slots (i.e., £ = integer- 7’); since symmetric pulses have been chosen for 
illustration in Fig. 1, their maxima are centered in the time slots for 
K odd, and lie on the time-slot boundaries for K even. Examine the 
(0, 7] time slot in Fig. 1 as typical; then with the above choice of 
pulse positions, the number of pulses contributing to ¢(¢) in each 
time slot equals K. Since each pulse can take on M different shapes, 
¢(t) can take on M¥* different shapes in each time slot; the same is 
true for v(t) of (41), thus demonstrating the representation of (43). 

It remains for us to express the pulse shapes r(f) and coefficients 
b; of (48) in terms of the signal pulses g(¢) and coefficients a, of (42). 
We give separate treatments for the cases K = 1 and K = 2, and 
extend these results to general K. 


4.1 Nonoverlapping pulses: K = 1 


The top portion of Fig. 1 shows digital phase modulation for which 
the signal pulses in different time slots never overlap; in this case, 


g(t) | = 0], (<0, t>T, (45) 
where 0 | represents the M-dimensional zero vector. Define 
[ein (t) ginlt) ... efoml)] Octs T. 
ad, = : (46) 
0, t<0, t>T. 


Then (41) and (42) may be written 


ut) = S aati — kD] (47) 


' 0] =,0,' is a vector of appropriate dimension (here M*) with all elements zero. 
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4T 





Fig. 1—Phase modulation for different signal pulse durations. Index k [eq. (42)] 
is shown near peak of each pulse. Also, for simplicity, same signal pulse is shown for 
each k. 7 = time slot duration or signaling period. KT = maximum signal pulse 
duration. Note different pulse center location for odd and even K. 


Comparing (47) and (43), the parameters of the latter are given as 
follows for nonoverlapping signal pulses: 


be= ay, rOl=aq@]; K=1. (48) 


4.2 Overlapping pulses: K = 2 


This case is illustrated in the second portion of Fig. 1. In the (0, 7] 
time slot, the k = 0, 1 pulses contribute. We have 


g(t) ] = 0], bad. tod. (49) 
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[ein (¢) eja(t) ... efom (t) | —T<tsT. 


O, poet ST 
Then 
of) = SL {awa — bT)]} faces: ae — (b+ )T)I) 


fax at — bT) J} X fauna b+ DT] 


= 2 bayX aers} tae — bT)| X ale — (e+ YT) I}, (1) 
where X denotes the (right) Kronecker product!?:!3 throughout.* 
Comparing the last line of (51) with (48), the parameters of the latter 
are given as follows when no more than two signal pulses overlap: 


be = a4,X ary rH] =q0]xXaqt—7T)]; K=2. (52) 


Since the elements of b; consist of all pairs of products of the elements 
of a; and az41, and since the a; are M-dimensional unit-basis vectors 


T Comparing (50) with (46), note that the definition of q(é) is different for different 
K; q(t) + 0 over the same interval in which g(t) may be nonzero. 

+ The second line of (51) follows from the first by the observation that the two 
scalars may be regarded as 1-by-1 matrices; consequently, their scalar product and 
their Kronecker product are identical. The third line follows from the second by the 
well-known result connecting ordinary matrix and Kronecker products as follows: 
Consider two arbitrary matrices A and B with elements a;;, b;;. Define the (right) 
Kronecker product by (Refs. 12 and 13): 


auB a2B ::: 
AXB= Q1B Ao2B eae 


The following results may be found in Refs. 12 or 13. 
For any matrices A, B, C, and D: 


AXBXC= (AXB)XC=AX (BX OQ), (A XB)’ = A’ XB’. 
For A and B the same size, and C and D the same size (possibly different from the 
size of A and B): 
(A+B)X(C+D)=AXC+AXD+BXC+BXD. 


Assume the matrices A1, Bi and Ao, Bz are dimensioned so that the ordinary matrix 
products A;-Bi and A2-Be exist (i.e., there are A columns of A; and rows of Bi, and 
» columns of A» and rows of Be). Then 


(Ai-Bi) X (Ae-Bo) = (Ar X Ae)- (Bi X& Bz); 
this result generalizes to 
(A1-B:-Ci) & (Aeg:B2-Ce) & (A3-B3-C3) 


= (Ai X A: X A;)+ (Bi X Bz X Bz)- (Ci X C2 X Cs) 
ec. 
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by (14) and (15), the b; are A7?-dimensional unit-basis vectors: that is, 
b; has one element unity and the remainder M? — 1 elements zero. 
Note from (52) and (50) that 


r(t)] =0], t<0, t>T. (53) 


We illustrate these relations for the binary case,!*4 in which only 
two signal waveforms gi(t) and ge(t) are transmitted. Then 


el=|%) |-[o]- fa FST. (54) 


ga(t) 0 
esas (t) 
eo) wees 
q(t) ] = (55) 
a t<—T, t>T. 


eila (t)ta1(t-1)] 
oilay (t)-+02(t-7)] 


eilos(t)tai(t—T)] | ? O<tsf. 
etlo2(t)t+92(t-T)] 
r(i)] = (56) 
0 
0 
ol» f=0, £> 7. 
0 
bz, = Lak? aP ] X Lats a®hi] 
Lat aah, af aft: a afl]. (57a) 
1 
[1 0] [0 1] 


(1000] | [0100] 
(0010) | fo001] 





The four possible waveshapes for ¢(¢) in the (0, 7] time slot are 
shown in Fig. 2. 
4.3 Overlapping pulses: General K 


The general case is treated by straightforward extension of the above 
method; Fig. 1 illustrates the phase modulation for K = 3, 4. The 
following expressions differ slightly for the even and odd cases; they 
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Fig. 2—#(t) forO0 < tS T, binary signaling, and K = 2. 


correctly reduce to the above results for K = 1, 2. 

















K-11 K+1 
ee oe ; 
ts 5 L.. tS 5 y i 
g(t) | = 0], 7 
K 
[e%a1(t) ejor(t) ... efom(t) | 
a A= eee ee Bee T; 
K K 
meaty eee eve . 
5 T<ts 9 f Bi 
até), = 
K-1 K+ 1 
9 7, t> 9 fi 
0, 
K Kk .. 
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K odd. 
(58) 


i even. 


K odd. 


K even. 
(59) 
K odd. 


K even. 


The parameters of (43) are: 


(K—1) /2 


(K—1) /2 . 
bi= Ux aus rO]= Tx alt—i2)]; Kodd. (60) 


t=— (K—1)/2 


K/2 K /2 
boe= IIx aay tOl= Hx aé¢-—77)]; Keven. (61)! 
i=—(K/2)+1 i=—(K/2)+1 
Note that 
r(t) ] = 0], (<0, t> fT. (62) 


V. PSK WITH NONOVERLAPPING, CORRELATED SIGNAL PULSES: K = 1 


Assume the signal pulses are confined to one time slot, as in (45): 


gt)}/=0], ¢50, t>T. (63) 


—- 


From (46) to (48) and (160) to (161), 
ef; (t) 


e192 (t) 


R(f)] = ie on ienft dt. (64) 


eso (t) 


We separate v(t) of (2) and (3) or (41) and (42) into line and continuous 
components: 


v(t) = e9%) = y(t) + 0, (t). (65) 


Using (48) and comparing (32) to (37) with (149) to (153), we have 
from (171) 


v(t) = ee > R( 7) ein ant /T. (66) 


— 2 


From (165), 
Pal) = palmROI? ES 5(f- FZ), 67) 


‘The symbol IIx used here and subsequently indicates a multiple Kronecker 
product: 


- . 
IIx A; = An X Anu X --- K Aw-i X& Ay. 
i=L 


The products in (60) and (61) contain K factors. 
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which may be written 


Pai(f) = qa ltsRa(f) + wolf) +--+ + waa)? 


2.8 #): 


Pi) = FRO; Se (Wy — whew) RD) 


From (162) and (164) 
We illustrate these general results with two examples. 


5.1 Independent, M-ary signal pulses 
Equation (32) now holds for all n ¥ 0; using (21), 


Wi 
W,=wiw= — -[w1We--: wu], n ¥ 0. 
was 
From (27), 
Wy 0 
Wo = Wa = ve 
0 WM 


Then only the n = 0 term remains in (69); 
] 1 
P..(f) = FRO) wa R*()] — 7 lw ROI. 
Writing out the matrix operations, (72) yields 


P,.(f) = 7 Cus | Ra(f) |? + we|Ro(f)|? +--+: + wm|Ra(f)\7] 


7 7 [wiki(f) + welte(f) + +++ + wu (f)|? 


1 M M 
sh aT zu 2 wi w;|R.(f) = R,(f) |?. 


(68) 


(69) 


(70) 


(71) 


(72) 


(73) 


This result has been given in Ref. 6. The line component is given by 


(66), its spectrum by (67) and (68). 
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5.2 Correlated, binary signal pulses: M= 2 


The matrix W, consists of the four upper left-hand elements of 
(23). The constraints of (19) and (20) or (25) and (26) render three of 
these functions dependent on the fourth. 


Wa= CWA — wi) | twee 
From (27) 
Wo(1, 1) = wi, (75) 
from (30) 
W,(1, 1) = W_,(1, 1), (76) 
and from (31) 
lim W,(1, 1) = wi. (77) 


n> 0o 


Equation (69) becomes 
P.(f) = plR() ~— Ral? Se PvMLWaCL, Y — wh}. (78) 


W,(1, 1) may be any discrete covariance function satisfying the con- 
straints of (75) to (77). The line component and its spectrum are again 
given by (66), (67), and (68) with two-component vectors, e.g., 


Piaf) = qelwik() + wiki?  o(f—-R)- 79) 


= 


The binary independent case is obtained by setting 





W,(1,1) = w, n 0. (80) 
Then (78) becomes 
P.(f) = “pr lRuf)— Ralf) |?. (81) 


This agrees with (73) for the binary case. 


Vi. PSK WITH INDEPENDENT, OVERLAPPING SIGNAL PULSES: K=2 


Assume that no more than two signal pulses overlap, as in (49): 


gij|}=0j], tS —T, t>T. (82) 
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R(f) is determined from (50) to (53) and (160) to (161): 


esl 9, (t)-+a,(t- T)} 


etl 91 (t)+ox (t—-T)} 
etl og (t) +9, (t-T)} 


ej{ 0g (t)+on (t-T)} 


T 
RNI= [ om sarst aon dt. (83) 


eit Ou—1 (t)-+oy, (t—-T)} 
etl on (t)+0;(t-7)) 


pilex O4ag0ST)) 
From (52), 
b, = ak, X an4i- (84) 
We determine the mean and covariance of b; from (32) to (87) and 


(149) to (153). We assume throughout this section that signal pulses 
in different time slots are independent, as in (70) and (71): 


Wo = wa; Wrz= wim, in| <0. (85) 
For the mean, since a; and a;41 are independent, 
8, = (bi) = WX Wy, 
= wil, (86) 


We introduce the notation that an integer exponent enclosed in square 

brackets denotes the Kronecker power,!? i.e., the Kronecker product 

of the matrix (or vector) with itself the indicated number of times. 
For the covariance (see footnote, p. 908): 


,(n) = (Di4n]- bx) 
= ((anin] X axing): (ai, x ars-1)) 
= ((arin] ax) X (an+n41]- aes). (87) 
Since 
®;(n) = ®,(—n), (88) 


we study only n = 0. We treat three cases. 
First assume |n| 2 2. All the a’s in (87) are independent; from the 
second line, 


@,(n) = (w] X w])-(w, X w) 


= wl?l-wPl, |n| 2 2. (89) 
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Next consider n = 0. From the third line of (87) and the fact that 
a, and ax%41 are independent, 


®,(0) = wi?! (90) 


Finally, for |n| = 1, from the third line of (87) we havet (see 
footnote, p. 908): 


(1) = ((anyi] & ax) X (ere) & anys) 
= (an, X (@rgi] X axpa) X ar+2]). (91) 


Since ax, az41, Aye are independent, 


®,(1) = w X wa X wi). (92) 
From (88), 
®,(—1) = w] X wa X wv. (93) 


From (86) to (89), (153) 1s satisfied, i.e., b, for widely separated 
time slots are uncorrelated.? Therefore, from (86), (165), and (171): 


nt) = qual XS R(g) fener (94) 


Pu(f) = qe lw -R(NIL? Daf — 22). (95) 


In comparing these results with (66) and (67), note that R(f) of (83) 
differs from R(f) of (64). 

Substituting (86), (89), (90), (92), and (93) into (162) and (164), 
the continuous spectrum is 


Pf) = GR), {(wH? — ww) 
sa (Wy, x Wad x w | — w IE] -wll)e rst 


+ (w] X wa Xw— wi?) wllet?rs7} R*(f) J. (96) 
We illustrate these results for the binary case. 
t The following vector relations may be established by inspection: 


a}-b,= a] Xb,=,b,X a]. 


*b.+n and bz also become independent as n — ~, because the b; are unit basis 
vectors (Ref. 10). 
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6.1 Independent, binary signal pulses: M=1 


R(f) in (83) now contains four components, the Fourier transforms 
of the components of (56). We have 


Ww 0 
w= [wi wi, we= | 2 (97) 
wl = [wi wiwe wie w3] (98) 
wi 0 0 0 
0 ww 0 0 
2] _ (We 
0 0 0 we 
Wt WiwWe Wiwe wiwe 
3 Dc 2,2 3 
W1iW2 WiWe WyiWe WiWe 
(2].wl2] = 
See, WiW. Wiwe wiws ww 0) 
wrws wwe ww ws 
we 0 wirwe 0 
2 2 
= WiWe2 0 WyW2 0 
Wire We a WL 0 wiwe OO wwe (101) 
0 wwe oO Ww 
wf] X wa X w= (Ww, X wa X w))’. (102)t 


Equations (98) to (102) substituted in (94) to (96) yield the final 
results for independent, binary signal pulses. 

To iulustrate, we write out a few terms of the matrix contained in 
{ $n (96), for the continuous spectrum in the binary case:!® 


{ da = wi(l — wi)(1 + 2w: cos 27 fT), 
{ }io = wiwret?27/? — wiwe(l + 2 cos 2rfT), 


Pee ae 


Writing out even the present simplest overlapping case produces an 
awful mess. Consequently, in the examples presented in this paper the 
digital computer is programmed to work directly with (96) or its 
generalization (116), performing matrix operations (both ordinary and 
Kronecker matrix multiplication) directly, rather than entering expres- 
sions such as (103). In this way, quite complicated cases involving 
multilevel signal pulses overlapping several time slots may be simply 
treated. 


(103) 


T See footnote, p. 908. 
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The general results for K = 2, (94) to (96), require the signal pulses 
to be less than two time slots in duration (82). Therefore, they must 
apply when the signal pulses are restricted to a single time slot, i.e., 
when the stronger condition (68) is satisfied, and must then reduce to 
the results for K = I, (66), (67), and (72). This reduction is demon- 
strated in Appendix C. 


Vil. PSK WITH INDEPENDENT OVERLAPPING SIGNAL PULSES: GENERAL K 


Finally, assume that the signal pulses occupy at most K time slots, 
so no more than K signal pulses overlap ; g(¢) is now given by (58). 
r(t) is given by (59) to (62), and R(f) by (160) and (161). 

The b; are given by (60) or (61); the mean and covariance of b, 
are found from (82) to (37) and (149) to (153), assuming throughout 
that signal pulses in different time slots are independent, as in (70) 
to (71). 

For the mean, 


B= (bx) = Ix ers) (104) 


by the independence of the a;. The limits over the [Ty, given in (60) 
or (61) for K odd or even respectively, extended over K factors in either 


case. Thus, 
6 = wi! (105) 


For the covariance (see footnote, p. 908): 


@,(n) = (Deyn ]- Dx) 
= (UI x @ninsi])- (II xaz+3)) 


= (IT x (Aitnti]axti)), (106) 


the [[, being taken over K factors as given by either (60) or (61). 
By stationarity, (106) is independent of k; consequently, for both even 
and odd cases 


(nm) = (burn by) 


~ (( Tx anss3)-( Tx as) ) 


= (Th (Anyi ]a) )° (107) 
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First, for n = 0, from the third line of (107) and the independence 
of the az, 


,(0) = Ix (az}.as) = whl, (108) 


Next consider n = 1. From the third line of (107) (see footnote, 
p. 915), 


®,(1) = (Tl (ai, X ais) 


= (ai) X | TI (ac]-ad)| X (ax+i]) 
= w,X wi] x wi], (109) 


the second line again following from the independence of the ax. 
For n = 2, proceeding as above: 


®, (2) = (ils (ai, X ais) ) 
= (ai) X (ix (aizi] X a)) X (ax+s]) 
= (a1) X (ix (ai, X ais+])) X (ax+2_]) 


= (as) x (as) X | Thx (acl-ad} X eal) X Caxsa) 


=wll X wi? & will, (110) 
By induction: 

®, (3) = wil & wok-3] x w |] 
2 ie 
0, (K — 1) = wit X wa X w | {4-1 ( ) 

®,(K) = wk] & w]I*), 
Next, from the second line of (107) and the independence of the ag, 
®,(n) = w]!Kl.wiKl) n= K. (112) 


This result is of course consistent with the final line of (111) (see foot- 
note, p. 915). Finally, 


®,(n) = @,(—7). (113) 
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From these results the line component is: 


vi(t) = wi*l. 3 R (3) eindrt/T. (114) 


Je AI 


Pal) = qolwe-R(NI? YS (t- 5) cas) 


The continuous spectrum is: 


P,(f) = 7 RU): (wart = Ewes) 


‘ 


7a 3 (wir x whk-n x w fir — w JA]. wIK])e—sn2rsT 


n= 


rary 


my 


+ > (wll & whE-" & wil — wt) wikletin2rsT} 
n=1 
‘R*(f) J. (116) 


If we set K = 2 in (114) to (116), these results agree with those of 
Section VI. If the signal pulses are restricted to a smaller number of 
time slots K < K, the present results reduce correctly to those ap- 
propriate for K; this reduction is shown for K = 1, K = 2 in Appendix 
C, but the general case is left as an exercise for the reader. 


Vill. EXAMPLES 


We now consider a number of examples of computing spectra of 
multiphase PSIX systems, subject to the following assumptions: 


(2) The number of phase levels is a power of 2, 
M = 2%, N an integer. (117) 
(iz) The M signaling pulses have a common shape, and the 
phase levels (peak phase modulations) are equally spaced. 
g(t) = ii [1 3--- (2M — 1)]9(2), (118) 


g(t) max = |. (119) 


Figure 3 shows vector diagrams of the peak phase modulations 
for M = 2, 4, and 8. We also assume that g(¢) is symmetric, 


g(t) = g(— t), K even, 


g(t + T/2) = g(—t + 7/2), K odd, (120) 
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m=8 


Fig. 3—Peak phase modulation for PSK with equally spaced levels. 


with maximum at 
t 


I 
© 
yy 
en 
< 
O 
— 


t= —, K odd. (121) 


(211) Signal pulses in different time slots are statistically indepen- 
dent, and all signal pulses are equally likely; 


l 1 
wm 11 = if 


where I is the identity matrix of order AZ. 


As a consequence of (118) and (122), we can show that P,(f) is 
symmetric, that is 


w]= Wa I, (122) 


P.(f) = P.(—/). (123) 


8.1 Rectangular nonoverlapping signal pulses 
If g(t) is a rectangular pulse (see Fig. 4) of duration n S T, 


T= a Bre ee 
g@)- 7 2 SSE 


0, otherwise. 





< 
O<a9sf, (124) 
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From (118) and (64), 
Raf) = lexpl i Fp @i— YD) — 1p AM eae 


= uae e-srfT y= 1,2,---,M. (125) 

The spectral density P,(f) can be computed from (68), (73), and 
(122) with R.(f) given by (125). The results can be shown to agree 
with those given in Refs. 4 and 6. Since this case can be treated as 
amplitude modulation and has been discussed extensively in Refs. 4 
and 6, we shall not discuss it further in this paper. 


8.2 Raised-cosine nonoverlapping signal pulses: K =1 


If g(t) is a raised-cosine pulse (see Fig. 5) of duration 7’, the pulse 
just fills the time slot, and 


] Qrt 

es — aes < 

51 cos |, 0O<isT 

g(t) = (126) 


0, otherwise. 


In this case, R(f)] may be evaluated either numerically or using 
Bessel function expansion (which again requires the use of a computer). 
For M = 2, 4, 8, and 16, we have evaluated as above 


P.(f) = P.(f) + P.(A), (127) 


with results shown in Fig. 6. We note that there is very little variation 
in either the discrete or the continuous spectrum when J is increased 
from 2 to 16. The tails of the spectrum are not shown in Fig. 6, although 
they are easily calculated down to the —60-dB level. 





Fig. 4—Square-wave signaling of duration yn S$ T. K = 1. 
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Fig. 5—Raised-cosine signaling with pulse duration T. K = 1. 


8.3 Raised-cosine overlapping signal pulses: K = 2 


If a raised-cosine signal pulse (see Fig. 7) just fills up two time slots, 


i 


g(t) = (128) 
0, otherwise, 


5[1 + cos % |, SS ROG Ss 


kK = 2, and two signal pulses overlap. 

In this case, P,,(f) and P,,(f) are given by (95), (96), (122). Using 
a digital computer, P,(f) has been determined for M = 2, 4, 8, and 
16, and the results plotted in Fig. 8. Again, it may be noted that there 
is not very much variation in either the discrete or the continuous 
spectra when the number of phase levels M is increased from 2 to 16. 

Comparing Figs. 6 and 8, we observe that the power in the line 
components in a PSK system with overlapping signal pulses is smaller 
than the power in the lines with nonoverlapping pulses. Also, for the 
same signaling rates, the principal portion of the continuous part of 
the spectrum with K = 2 is narrower than that with K = 1. The tails 
of P,(f) were calculated down to —60 dB even though they are not 
shown in Fig. 8. 

In Fig. 9, for M = 4, and raised-cosine wave signaling, we plot the 
spectral density P,(f) when the amount of overlap is increased from 
zero to one time slot. Specifically, we assume that 


1 at 
= Ride = < <g< 
5 [1+ 008 | bf <1 s 67, 10,9 = 6-=- 1, 
g(t) = (129) 
0, otherwise. 


Note that K = 1 when 6 = 0.5, and K = 2 when 0.5< 6S 1. 
Figure 9 shows that there is a gradual reduction of power in the line 
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ve POWER IN LINE COMPONENTS 


3.63 x 10-1 | 3.30x 10-1 | 3.23 x 10-1 
6.60 x 10-2 | 6.99 x 10-2 | 7.07 x 10-2 | 7.09 x 10-2 
2.68 x 10-3 | 1.47x 1073 | 1.24x 1073 | 1.18x 1073 


471x 10-5 | 2.15x 10-4] 2.63x 10-4 | 2.75 x 10-4 


SPECTRAL DENSITY P, (f) /T 





0 1 2 3 
NORMALIZED FREQUENCY fT 


Fig. 6—Spectral density of binary, quaternary, octonary, and 16-ary PSK systems 
with raised-cosine signaling and pulse duration 7. K = 1. Although the values of the 
spectral density for M = 2, 4, 8, and 16 are shown to be the same over certain seg- 
ments of the range of f, they are usually different from each other, but this difference 
is not large enough to be shown on the linear scale in Fig. 6. 


components when the overlap is increased from zero to one time slot. 

There is also the narrowing of the principal portion of the spectrum. 
Decrease of the carrier component [7 = 0 term in (68) or (95) | when 

6 1s Increased may result from the fact that the phasors in Fig. 3 spend 
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Fig. 7—Raised-cosine signaling with pulse duration 27’. K = 2. 


increasingly less time in the neighborhood of 0°. There are also smoother 
transitions between the phasors with increasing 8; this probably ex- 
plains the decrease in width of the main part of continuous spectrum. 
The narrowing of the spectral density and the reduction of the discrete 
spectral lines may or may not indicate better interchannel and inter- 
symbol performance, but this remains a subject for future investigation. 


IX. SUMMARY AND CONCLUSIONS 


Matrix methods have been used to express the spectral density of a 
sinusoidal carrier phase modulated by a pulse train with a finite pulse 
duration. Arbitrary pulse shapes may be used for M-ary digital signal- 
ing and they may overlap over a finite number of signaling intervals. 

If the pulse duration is one signaling interval, our results give the 
spectral density even though successive symbols are correlated. If 
the pulse duration is more than one signaling interval, we assume that 
symbols transmitted during different time slots are statistically inde- 
pendent; the case of correlated symbols may be considered by ex- 
tension of the present work, but the required statistical description of 
the modulation a; will be more complicated. For example, if K = 2, 
in addition to (ax |), (ax] X a1), we need to know (a;] X az] X a,X an), 
the fourth order statistics of a,, to determine the spectral density from 
our methods. 

The spectral density has a compact Hermitian form suitable for 
numerical computation by a digital computer. The associated Hermi- 
tian matrix is a function of only the symbol probabilities, and the 
column matrix associated with the Hermitian form is the Fourier 
transform of certain time functions related to the signaling pulses. The 
computations presented in this paper for binary, quaternary, octonary, 
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1.2 POWER IN LINE COMPONENTS 


NUMBER OF PHASE LEVELS 


1.0% 1.13 x 10-2 


1.03 x 10-4 | 1.12 x 10-4 | 1.14 x 1074 


1.88 x 10-7 | 2.58x 10-7] 2.78x 1077 


0.8 


SPECTRAL DENSITY P,(f) /T 
(o) 
ro?) 


0.4 


0.2 





NORMALIZED FREQUENCY fT 


Fig. 8—Spectral density of binary, quaternary, octonary, and 16-ary PSK systems 
with raised-cosine signaling and pulse duration 27. K = 2. Although the values of 
the spectral density for WZ = 2, 4, 8, and 16 are shown to be the same over certain 
segments of the range of f, they are usually different from each other, but this dif- 
ference is not large enough to be shown on the linear scale in Fig. 8. 


DIGITAL PM SPECTRA 925 


1.2 POWER IN LINE COMPONENTS 


OVERLAP PARAMETER 
0.500 0.571 0.667 0.800 | 1.000 | 


3.30 x 10-1 
















2.64 x 10-' | 1.88 x 10-1 | 1.05 x 10-1 | 3.73 x 10-2 











6.99 x 10-2 





6.69 x 10-2 | 5.52x 10-2 | 3.33 x 10-2 | 1.12 x 10-2 















1.47 x 10-3 | 5.32 x 10-5 | 2.95 x 10-4 


1.30 x 10-6 


4.92 x 10-4 
3.57 x 10-6 


1.03 x 10-4 




















2.15 x 10-4 | 9.45x 10-5 1.88 x 1077 





0.8 


0.6 


SPECTRAL DENSITY P, (f) /T 


0.4 


0.2 


NORMALIZED FREQUENCY fT 


Fig. 9—Spectral density of quaternary PSK system with raised-cosine signaling 
and pulse duration 287, 0.56 S 6 S$ 1. When B = 0.5, K = 1, the pulse fills up just 
one time slot and the spectral density is as shown in Fig. 6. When 6 = 1, K = 2, 
the pulse fills up just two time slots, and the spectral density is as shown in Fig. 8. 
Although the values of the spectral density for 8 = 0.5, 0.571, 0.667, 0.8, and 1.0 
are shown to be the same over certain segments of the range of f, they are usually 
different from each other, but this difference is not large enough to be shown on the 
linear scale in Fig. 9. Note that when 6 = 0.5, 0.571, 0.667, 0.8, and 1.0, 1/8 = 2, 
1.75, 1.5, 1.25, and 1.0. 
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and 16-ary PSK systems do not consume very much computer time 
and are quite inexpensive. 

Extraction of a timing wave is often essential in the detection and 
regeneration of digital signals. In a self-timed digital system not 
containing a timing wave, the wave is extracted from an incoming 
pulse stream by a nonlinear processing of the signal. For example, the 
incoming pulse stream may be passed through a square-law rectifier 
and a linear narrow-band filter to get the timing wave. We would like 
to note here that the methods given in this paper can be extended to 
determine the spectral density of the output of a nonlinear device 
whose input is a random time pulse train of the form (16) or (42). The 
results presented here apply to the particular nonlinearity exp j(-); 
other nonlinear functions are treated in a similar manner. 
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APPENDIX A 
Spectra of Complex and Real PSK Waves 


From (1) written in complex form, the covariance of the real wave 
a(t) is 


b,(r) = &.(t, 7) 
= Le2rferh, (7) tea ferpi(z) + ei rloreitrseth,,«(t, 7) 
ferdrtereitrse@en(t, 7)}, (180) 
where the cross-variance is given by 
P,y*(t, rT) = (u(t + r)v(t)) = (etletnte), (131) 


Now ®,,+(t, 7), with 7 regarded as a parameter, is periodic in ¢ with 
period 7; 3 
b,.*(t + T, 7) = Byp*(t, 7). (132) 


This follows from the assumed strict stationarity of the s, of (2). Then 
we may write 


,,+*(¢, 7) = 3 On(remertlt, (133) 


nr=— 
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Substituting in the time average quantity of (130) and interchanging 
the order of summation and integration, 


ein fet@,,*(t, 7) 


. A 
> On{T) lim = es2rl2fet(n/T)] tt 
Tes Av o@ 2A Pay 


S y,(r) lim sin 2n[2f- + (n/T)]A 


2, o7) jim “Spar (8M 


The im — 0 if 2f. + (n/T)# 0 for every integer n, 1.e., if twice the 
Aw” 


carrier frequency is not a precise multiple of the modulation fre- 
quency 1/7. Under this condition, (180) and (134) yield (5). 
Next, assume that P,(f) is strictly band-limited; 


Pf) =0, [fl 2 fe. (135) 


Then P,(f — f.) and P,(—f — f.), the two terms appearing in (5), 
do not overlap. Now define 


v(t) = eirFety(t), (136) 
Then 
B(t, 7) = e?*Fe7h, (t, 7); Bir) = et Fe7h, (7). (137) 
Pf) = Pus — fs Prof) = Po(-F— fo). 
Therefore, subject to (135) 
f = 0. 
Pp, (f) aos 0, 
J 2 23 
(138) 
P,+(f) = 0, 
f2%: 
Now the cross-variance of y(t) and »*(t) is 
belt) = Gall, 8) = Ole THT 7), (139) 


where ®,,* is given by (181). The Fourier transform of (139) is the 
corresponding cross-spectrum P,,*(f). Since the two spectra P,(f) and 
P,*(f) do not overlap by (138), the cross-spectrum P,,*(f) must be 
identically zero,!* and hence also the cross-variance; 


6,*(7) = 0, P.(f) =0 for |f| = fe. (140) 


928 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1974 


Substituting (139) and (140) into (130) and taking the Fourier trans- 
form, 


Pf) =0 for |f| 2 f. (141) 


While P,(f) will never be strictly zero, it will eventually fall off so 
rapidly with increasing | f| that the spectral relation of (5) will be a 
good approximation when f, is large enough so that the two terms do 
not overlap appreciably, even if 2f, is an integral multiple of 1/7. 

Consider now the exceptional case of a low carrier frequency that is 
an integral multiple of half the modulation frequency: 


of. + 7 = 0. (142) 


Then in (1384) the lim — 1 for the no term, and —0O for all other 
A-© 
terms as before. Therefore, (180) and (134) yield 


@2(7) = cePrFe"T Dh, (7) + g-azer(7) | 
+ Le-i2r Se Bir) + yor (r)], 2f.7' = integer. (143) 


The exceptional case requires the evaluation of the additional quantity 
y-oy.r(7), 2f-T' = integer, where ¢ is defined in (133). This may be 
done similarly to the evaluation of ®,(7) performed in the text, but 
it is not undertaken here. 

Finally, consider the wave of (9), with #(t) given by (2), with the 
carrier and modulation phases ¢o and ¢) independent random variables, 
independent of the modulation. Equation (130) is replaced by 
Bi(r) = Zler"Forh, (7) + err her@y (2 

+ e727 fe7(ei2bo) (gifts eit feth, .«(E, T) 
feriteler(e-itds)(enitrseta)eHTBea(t, )], (144) 
where v, ®,, and @,,* remain as given in (3), (7), (8), and (131) for 
to = 0. The last two lines of (144) vanish if any of the three quantities 


eHFIB.-(t, 7), (ei), (eitr/et) (145) 
vanish. Therefore, 


P.(f) = 2P.(f — f.) + ¢P(—-f — fi) (146) 
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if any of the following conditions are satisfied: 


2f.T A integer. 
2f.T' = integer, ¢ois uniformly distributed over an interval 


w, 27, 37, °°. (147) 
2f.T' = integer, tf) is uniformly distributed over an interval 
7 T 9 T 
va ie tae 250° ye ed head 


The first condition is, of course, that of (5). The last two conditions 
show that suitably randomizing the phase of either the carrier or the 
modulation suffices to yield the simple spectral result of (5) when 
2f,.7’ = integer. 


APPENDIX B 


Spectrum of a Baseband Pulse Train With Different Pulse Shapes 


We determine the spectral density of (43), 
v(t) = Yo byt(t — kT), (148) 
k=—0o 


by the vector analog of the corresponding derivation for a train of 
equally spaced pulses with similar shapes. In this appendix, the wave- 
forms r for different time slots (i.e., different k) may overlap, although 
they do not overlap in the applications of this paper. Also, b; may be 
complex in this appendix. It is real in the main part of the paper. 

Let b; in (48) and (148) be wide-sense stationary. Then, using the 
notation of (35) to (40): 


G = (bx). (149) 
®,(k — 1) = (b;]-bi). (150) 
Pi(f) = Ye Ps(n), (151) 
By(n) =f ByNermmnas, (152) 
Further, assume that b;,, and b, become uncorrelated as n —> «: 
lim @,(n) = (0) = 6]-B*. | (153) 
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From (148), 





®,(t, 7) = (u(t + 7)v*(t)) 
= > = r(t + r — kT) -@,(b — 1)-r*(t — UP)). (154) 
Since 
&,(t + 7,7) = &,(t, 7), (155) 
we have 
ee eae eee Ut (156) 
v v\Yy T rye vl, . 


Substituting (154) in (156) and making the substitutions k — 1 =n 
andi —IT =, 





1 2 00 —(l—$)T ~ 
Spee ee i; r(t + r — nT) -@,(n)-r*(t) Jat 
DP yee are a 
o 7. iz r(t + + — nT) - y(n) -1*(t) Jd. (157) 





The spectral density of (148) is now the Fourier transform of (157): 


P,(f) = iC b, (re? dr 


= ; y | ; | © ease r(t + r= nT) 


ae -,(n)-r*(t) Jdtdr. (158) 


The treatment differs slightly from that for the scalar case, because 
the order of the factors in the matrix products may not be changed. 
Setting ¢’ = ¢ + 7, (158) becomes 


P.(f) = = = Hf en serft! r(f! — nya -@,(n) 


n= 


| i * gtiteft y¥(t) Jd}. (159) 
Define 
RN] =R(N/ = fo ear Tat (160) 


i.e., the elements of r(¢) and R(f) are the pulse shapes and their 
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Fourier transforms, respectively. In applications throughout the re- 
mainder of this paper, 


ri)|]=0], ¢S$0, ¢t> 7, (161) 


and the integral in (160) therefore takes the limits /o7. However, the 
restriction of (161) is not necessary in this appendix. Substituting 
(160) into (159) and using (151), 


Pf) = FRU) Po(s7) RAI (162) 


our desired result. Eq. (162) reduces correctly to the scalar case.!” 

It is convenient to separate out the line component of v(t). As in 
(39) and (40), the line and continuous components of (151) are, using 
(153), 


Pu(f) = B1.gF Daf — 2). (163) 
P.(f) = De #"(Ho(n) — 81-8"). (164) 


The line and continuous spectral components of (148) are found by 
substituting (163) and (164) into (162). For the line component 


P.(f) PRD Og RIX 3(t- 5) 


| 


mIGRMI ¥ o(t- 5). (165) 


The line component is composed of de and sine waves at the signaling 
rate and its harmonics. 
Finally, taking the expected value of (148) and using (149), 


(v(t)) = 8: pas r(t — kT)). (166) 
Then 
@(,(t, r) = (v(t + 7)){v*(t)) 


2 2 r(t + 7 — kT)-6]-BF-r*@ — LT) ]. (167) 


I 
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Comparing with (154) and proceeding as above, 





B(,)(7) = a. i r(t-+ + — nT)-B)-G*-r*()Jdt. (168) 


Consequently, (168) is a periodic function of 7 with period 7. Compar- 
ing (157) and (168) as |r| — ©, using (153), and assuming that r(t) 
eventually falls off for large |¢| [the stronger assumption (161) applies 
throughout the rest of this paper ], 


®,(7T) > ®iy)(7) as [rl] mo. (169) 
Therefore, we may write (148) as 


v(t) = vi(t) + v,(t), (170) 
where, from (166), 


vi(t) = B- Oy r(t — kT) J 


= 78: ees R ( 7) len”, (171) 
and the spectral densities of v:(t) and v,(t) are given by (165) and by 
substituting (164) into (162), respectively. Equation (171) contains 
phase information lost in (165). 

The assumption of (153) follows if the vector coefficients b; become 
independent for widely separated time slots. For the applications of 
this paper, the b, are unit basis vectors, i.e., each b;, has a single 
component equal to unity and all other components zero. In this 
special case, the assumption of (153) that b; in widely separated time 
slots are uncorrelated also guarantees independence.” 


APPENDIX C 
Simplification of PSK Spectra for K = 2 for Nonoverlapping Signal Pulses 


We demonstrate that the K = 2 results, (94) to (96), which apply 
subject to condition (82), specialize to the K = 1 results (66), (67), 
and (72) when the stronger condition (63) is satisfied. 

We distinguish the different R’s in Sections V and VI by appending 
subscripts K, denoting (64) and (83) as R; and Ro, respectively. Then 
if, in Section VI, (82) is replaced by the stronger condition (68), (83) 
becomes 


R(f)] = Ri(f)] X 1], (172) 
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where R; is given by (64) and 1 is the M-dimensional vector with all 
components unity defined in (24). 

Considering first (94) and (95), using (26) and (172) (see footnote, 
p. 908): 

wi?l-Re(f)] = (w, X w)-(Ri(f)] X 1) 
= (w,Ri(f) ]) X (w1]) = w Riff). (173) 

Substituting (173) into (94) and (95) yields (66) and (67). 

Next, substituting (172) into (96), we evaluate the following resulting 
products using (26) and (28) (see footnote, p. 908): 


R2(f),- wi"? -Ro(f) ] 
= (Ri(f), X 1)-(wa X wa): (Ri(f)] X 1) 
= (Ri(f)-we-Ri(f)]) X @; wa-1)) 
= Riff) wa-Ri(f)]. (174) 


Ro(f)w J? -wi?].Ro(f) ] 
= (Ri(f), X 1)-(w] X w])-(w, X mw): (Ri(f)] X 1) 
= (Ri(f)-wJ-w-Ri(f)]) ¥ (y wl-w,1)) 
= (Ri(f),-w])(w,-Ri(f)) 
= |w,Ri(f)]|?. (175) 


R.(f)-(w, X wa X w]) RAS) 
= (1 X Rif), X 1)-(w X wa X W))- (RUT X 1 XD) 
= (1-w-RUN) X Ri(fwael) X ew)-D 
= (w;Ri(/)) (Ril), w) 
= |wRA(/)II?, (176) 


Substituting (174) to (176) into (96) yields (72). 
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Consecutive Delta Modulations of 
a Speech Signal 
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We consider a communication link in which a band-limited speech 
signal is delta-modulated, detected, and filtered by a low-pass filter, and the 
analog output ts delta-modulated again with an tdentical encoder. We are 
concerned with the correlation C' between equal-length bit sequences, desig- 
nated {b} and {B}, that result from the two stages of delta modulation. We 
study C as a function of the sequence length W; the starting sample T 
in {b}; the time shift L between {b} and {B}; the signal-sampling fre- 
quency Ff; and a parameter P(=1) which specifies the speed of step-size 
adaptations in the delta modulators. (P = 1 provides nonadaptive, or 
linear, delta modulation.) 

Computer stimulations have confirmed that for small time shifts L and 
for statistically adequate window lengths W, C'1s a strong positive number 
(0.4, for example). Moreover, the C function tends to exhibit a maximum 
Cmax at a small nonzero value of L (between 1 and 5, say) reflecting a delay 
entroduced by the low-pass filter preceding the second delta modulator; 
and when W 1s on the order of 100 or more, the dependence of Cmax on the 
starting sample T is surprisingly weak. Also, in the range of F and P 
values included in our simulation, Cmax increased with F and decreased 
with P. Finally, the positive C values for small L are retained even when 
the delta modulators are out of synchronization in amplitude level and 
step size, as long as the delta modulators incorporate leaky integrators and 
finite, nonzero values for maximum and minimum step size. 

With a gwen T, the C(L) function can exhibit significant nonzero 
values even for large L. However, these values are both positive and nega- 
tive; and if correlations are averaged over several values of T', the average 
C(L) function tends to be essentially zero for sufficiently large L (L = 100, 
say), while still preserving the strong positive peaks at a predictable small 
value of L. This observation is the basis of an interesting application 

937 


where the value of C 1s used to determine whether or not two digital codes, 
appearing at different points in a speech communication system, carry 
identical speech information. 


1. THE PROBLEM — 


Consider a speech signal subjected to two successive stages of delta 
modulation, with an intermediate stage of low-pass filtering, as shown 
in Fig. 1. A previous paper! has studied how signal quality degrades 
as a function of the number of delta modulations. The present paper 
is concerned with the amount of correlation that exists between the 
bit sequences {b} and {B} from the two (identical) delta modulators. 
Specifically, we employ computer simulations to study the correlation 


1 TiW 
Cc = Ww & O:Ber. 


It is assumed that {b} and {B} are zero-mean sequences with equi- 
probable +1 entries. Apart from being a function of the window 
duration W and time shift L, the correlation C’ will also depend on the 
signal-sampling frequency F and a parameter P specifying the step- 
size logic used in the delta modulators. The delta-modulator simu- 
lations are described in Section IT and the properties of C' are described 
in succeeding sections. 

The studies reported in this paper were prompted by an interesting 
potential application where the value of the correlation C’ would be 
used to determine whether or not two digital codes (appearing at 
different points in a speech communication network) carry the same 
speech information. More specifically, we were considering a telephonic 
system that incorporated digital and analog signal terminals capable 
of being interconnected via a common switching network. The problem 
was to determine whether digital terminals communicating with each 
other (in other words, handling the same speech information) could be 
detected by digitally correlating the signals of each digital terminal 
with the signals at other digital terminals in the system.?? The digital 
coding under consideration was delta modulation, and the results of 
this paper indeed suggest that the detection of communicating termi- 
nals should be possible on the basis of appropriate bit correlations. 


BAND — LIMITED BIT BIT 
SPEECH INPUT DELTA SEQUENCE eae DELTA SEQUENCE 
R 
% MODULATOR {o} ean (eer MODULATO {ey 


Fig. 1—Block diagram of the simulated speech communication system. 
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Fig. 2—Schematic diagram of an adaptive delta modulator. 


li. SIMULATION DETAILS 


The delta modulator utilized in our simulations is schematized in 
Fig. 2 and is the same instantaneously adaptive delta modulator 
(ADM) discussed in Ref. 4. Basically, it is described by the equations 


b, = sen (X, — Y,-4), 
Y,= | Goes a A,*0;, 
and 
A, = A, Por ort, 


where X, is the amplitude of the input sample r, and Y,_; is the ampli- 
tude of the latest staircase approximation to it. The parameter P 
(21) automatically increases step size when Y is not tracking X 
fast enough (6, = b,-1), and decreases it when Y is hunting around 
X (b, = — 6,1). Nonadaptive or linear delta modulation (LDM) 
corresponds to the special case of P = 1. 

The speech signal is a 1.5-second male utterance of “‘Have you seen 
Bill?” that is band-limited to 3.38 kHz. The sampling rate, unless 
otherwise noted, is 60 kHz. A plot of the speech waveform appears in 
Fig. 3, where a number at the right of a line represents the last 60-kHz 
sample in that line. The original signal samples are quantized to a 12-bit 
accuracy, and have integer amplitudes in the range —2!! to +2!!. 
Finally, the low-pass filter is a programmed recursive filter with an 
18-dB/octave roll-off. This seems to represent adequate filtering for 
toll-quality speech reproduction using ADM at 60 kHz. 


ll]. DEPENDENCE OF CORRELATION ON TIME SHIFT 


Figure 4 shows the dependence of C' on the time-shift L for two 
different values of starting sample 7’. It is interesting to observe that 
both the functions show a maximum at L = Lmax = 4. Even more 
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Fig. 3—The speech waveform of ‘‘Have you seen Bill?” 


interesting is the fact that respective values of CXZ(L), the correlation 
between the speech input X and the low-pass filter output Z, are also 
maximized (as determined in a separate simulation) at L = 4. It 
would seem that the nonzero value of Lmax 1n Fig. 4 is to be attributed 
to the delay introduced by the low-pass filter. Actual values of Cmax 
and Lmax depend on the short-term speech spectrum and the nature of 
low-pass filtering, as determined by the parameters 7’ and F (see 
Tables I and II). It is a general result, however, that the C(L) function 
always shows a unique, strongly positive maximum value at a small 
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Fig. 4—Dependence of correlation C on time shift D. 


value of L. Secondary peaks at large values of L tend to be less unique, 
and they tend to be randomly positive and negative depending on the 
part of the speech utterance being considered, as determined by T. 


IV. DEPENDENCE OF MAXIMUM CORRELATION ON SAMPLING FREQUENCY 


Table I indicates a tendency for Cmax to decrease with decreasing 
sampling frequency. This may be ascribed to the fact that, at a lower 
sampling rate, delta modulation provides a cruder approximation to 
the input signal. The bits, therefore, carry more signal-independent 
noise information, and they have corresponding random properties 
that cause a decorrelation between {6b} and {B}. 


Table |—- Dependence of maximum correlation Cyax 
on sampling frequency F (P = 2, W = 1000) 
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Table Il1— Dependence of correlation C on starting sample T 
and time shift L (W = 150, P= 2, F = 60 kHz; 
numbers in parentheses are values of Lax) 


- 0. Tee 10 100 1000 —- 10000 

7425 (4) 0.27 0.69 0.20 0.08 0.04 0.01 
17425 (4) 035 0.48 024 -009 —-004 —Ol1l 
27425 (5) 0:23 0.47 011 -003  —0.03 003 
37425 (4) 015 037 4-011 —0.16 —0.08  —0.08 
47425 (2) 0.29 0:33 031 —013 —O11 —0.08 
57425 (4) 037 0.43 011 0.04 0.21 0:13 
67425 (1) 039 0.44 0:37 0.27 0:22 —003 

Average of C 

values (over T)| 0.29 ~—0.46 0.18  —0.01 0.03 0.02 


V. DEPENDENCE OF MAXIMUM CORRELATION ON STEP-SIZE MULTIPLIER P 


Table III demonstrates how Cmax tends to decrease with increasing 
P. Larger values of P increase the high-frequency excursions of the 
staircase function Y. These are filtered out by the low-pass filter. 
This leads to lesser correlation between the filter output Z and the bit 
sequence {b} and, thence, to a decorrelation of {B} and {b}. 


VI. DEPENDENCE OF MAXIMUM CORRELATION ON WINDOW LENGTH 


Our results so far have tacitly assumed window length values that 
represent bit sequences whose durations are of the order of a few 
milliseconds. Figure 5 shows C' explicitly asa function of W. It is seen 
that very stable indications result with W in the order of 1000, although 
values close to a respective asymptote are sometimes reached for W 
values in the order of 100. In fact, a window length of W = 10 is seen 
to be sufficient, for all values of T in Fig. 5, to bring out the strong 
positive nature of Cmax. The convergence of the three curves in Fig. 5 


Table Il!— Dependence of maximum correlation Cy. on 
step-size multiplier P (F = 60 kHz, W = 1000) 


T 37425 37000 
Va Dmax Cmax Dmax Cmax 
1.0 4 0.91 4 0.89 
1.5 4 0.66 4 0.62 
2.0 4 0.34 4 0.44 
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Fig. 5—Dependence of correlation C on window length W. 


is not at all surprising. Note that, by definition, C should indeed be 
independent of 7 for W > ~. The results of Fig. 5 were based on a 
search for Cmax In the range 0 S L S 10. Except for W = 1, unique 
maxima were noted at nonzero values of L. For W = 1, the value of 
Cmax WaS surprisingly constant in the range 0 S L S 10, the constant 
value being +1 for one value of 7, and —1 for the other two. 


Vil. DEPENDENCE OF MAXIMUM CORRELATION ON WINDOW LOCATION 


As seen in the last section, Cmax iS a significant function of the start- 
ing sample for finite W. Table II shows the values of Cmax for seven 
equally spaced values of 7. The average value of Cmax iS 0.46 and the 
standard deviation is only 0.10. Notealso that C values for large ZL 
are smaller in general, and the effect is more noticeable when corre- 
lations are averaged over 7’. This is because the positive Cmax values 
always add up, while C' values for large L, being randomly positive or 
negative, tend to average out to values close to zero. 

At least one interesting application of the preceding observations has 
been suggested.?:? Suppose the second delta modulator has several 
potential speech inputs including the input Z resulting from X. The 
function C would then assume the strong positive values of Table II 
only when the input to the second delta modulator is indeed the DM 
version of the speech X ; and it would show values of C — 0 if the input 
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Table I1V—The effect of unsynchronized delta modulators 
(T = 37425, P = 2, F = 60 kHz, W = 1000. Values in 
parentheses are for W = 150) 


Initial Conditions 
> => Integrators 


Limits 
Cas Yi Y, Ai As Step Size Lines C max 
I 0 01 1 Perfect (0, 0) 4 (4) | 0.84 (0.37) 
II 0 0 1 1 Leaky (25, 250) 2 (8) | 0.48 (0.83) 
Itt 0 -—-50 1 —10 Leaky (25, 250) 2 (3) | 0.47 (0.76) 
Iv. | 0 —50 1 —10 Perfect (0, «) 5 (1) | 0.11 (0.16) 


was an entirely different speech signal* (possibly due to a different 
speaker). This effect will be more pronounced if the averaging process 
indicated in Table II is carried out. We are suggesting, in other words, 
a means of determining whether or not two digital DM codes, appearing 
at different points in a speech communication network, carry the same 
segment of speech information. The basic recipe is a DM bit correlator 
with a window of 0.1 to 1 ms, and a window location T that seems to 
be quite uncritical, especially when time diversity (averaging over 7’) 
is possible. 


Vill. EFFECT OF UNSYNCHRONIZED DELTA MODULATORS 


In practice, the two delta modulators in Fig. 1 can be unsynchronized 
in amplitude Y and step size A when either or both of them are in some 
kind of a transient state of operation. It is an interesting result of our 
study that the strong positive values of Cmax are retained even during 
such asynchronous periods, provided the delta modulators operate with 
a leaky integrator, and with finite and nonzero limits on step size. 
Leaky integration decreases the effect of amplitude history and, hence, 
the effect of amplitude asynchrony. Finite and nonzero limits on step 
size provide potential meeting points for the two step-size sequences, 
although they may begin with a different starting value. 

In Table IV, Y; and Y> represent initial amplitudes for the two delta 
modulators, while Ai and A, are the initial (signed) step sizes. The 
step-size limits, 250 (maximum) and 25 (minimum), include a signifi- 
cant range of step sizes that are called for in the adaptive delta modu- 
lation of speech (with F = 60 kHz, and with signal amplitudes in the 
range —211 to +211).4 Finally, the leaky integrators of Table IV 
leaked 5 percent of signal amplitude in a sampling period. 


"This situation is hypothesized to be equivalent to the case of large L. 
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Fig. 6—Dependence of correlation C on time shift 2 = —example of unsynchronized 
delta modulators. 


Note that Table IV shows that leaky integration and finite, non- 
zero step-size limits are imperative in the asynchronous case (rows 
III and IV, Table IV) to preserve a strong positive Cmax; they are 
also desirable to boost the value of Cmax in the synchronous case 
(rows I and II, Table IV). (The boost is quite significant for W = 150). 
A separate simulation showed that finite (nonzero) step-size limits 
and leaky integrators were effective only when employed in unison; 
and in one study of C as a function of L, they also sharpened the peak 
at Dmax (see Fig. 6). 

Finally, Table V is a counterpart of Table II for the case of un- 
synchronized encoders. The step-size limits are 25 and 250, the leak 
is 1 percent in a sample duration, and P = 1.5. (The last two numbers 
are probably more representative than the corresponding values in 


Table V— Dependence of correlation C on starting time 7 and 
time shift L with unsynchronized delta modulators 
(W= 10, P=1.5, F = 60 kHz) 


L 


i 0 Dixy 10 100 1000 10000 
7425 —0.4 0.0 0.0 0.0 0.0 0.0 
27425 —0.4 0.4 —0.4 —0.4 —0,.2 0.2 
47425 0.4 0.4 0.4 ~—0.4 —0.2 —0.2 
67425 0.6 0.6 0.6 0.6 0.6 0.0 
Average of C 
values (over T) 0.05 0.35 0.15 —0.05 0.05 0.00 
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Table IV.) Finally, we have reduced the window duration to W = 10. 
This results in obviously crude C(Z) functions for a given beginning 
sample 7’. But, as in Table II, when C(L) values are averaged over 7, 
the resulting C function shows a clear tendency to decay to near-zero 
values for LZ 2 100. The values of Cmax in Table V represent largest 
values as seen in a finite search (0 S L S 5). None of these was a 
unique maximum, which is possibly due to the insufficient duration 
(0.16 ms) of the window, W = 10. 
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Optical Waveguides With Very Low Losses 


By W. G. FRENCH, J. B. MACCHESNEY, 
P. B. O’;CONNOR, and G. W. TASKER 


(Manuscript received April 25, 1974) 


Low-loss optical fibers may be necessary for economical optical 
transmission systems. We have developed fibers that exhibit losses 
of less than 2 dB/km, at 1.06 um. The fibers were made by a chemical 
vapor deposition (CVD) technique that employs simultaneous re- 
action and fusion to a clear glassy core material. 

Two fiber compositions have been used. In the first fiber, a GeOo- 
doped fused-silica core is deposited inside a fused-quartz tube that 
acts as the cladding after the tube is collapsed into a rod and pulled 
into a fiber. 

Figure 1 shows the loss spectrum of a fiber made in this manner. 
The fiber is 723 m long and has a core approximately 35 wm in diameter. 
The numerical aperture is 0.235. The loss decreases by approximately 
\~4, the expected Rayleigh scattering dependence, to a minimum Just 
under 2 dB/km at 1.06 wm. Hydroxyl-ion-related absorptions at 0.72, 
0.88, and 0.95 um are low, amounting to less than 10 dB/km at 0.95 um. 
We believe that the OH impurities causing these absorptions are due 
to siloxane present in the SiCl, starting material. This can be removed 
by fractional distillation, and loss peaks due to the hydroxyl-ion- 
related absorptions as low as 2 dB/km above background at 0.95 um 
have been observed in similar fibers. This process has been used to 
produce GeO:-doped fibers with numerical apertures as high as 0.35, 
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Fig. 1—Loss spectrum of a fiber waveguide with a GeO2SiOe core and a SiOz 
cladding. 


and lengths up to 1.2 km. The length is presently limited by the avail- 
able fiber-drawing facilities. 

Figure 2 illustrates the loss spectrum of a second type of fiber 
consisting of a pure fused-silica core and borosilicate cladding. In this 
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Fig. 2—Loss spectrum of a fiber waveguide with a SiOe core and B2O;-SiO2 
cladding. 
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3—Loss spectrum of a fiber waveguide with a graded-refractive-index core 
BO 00 and borosilicate cladding. 


case, both core and cladding were formed by CVD with simultaneous 
fusion inside a fused-quartz tube, which will be described in an article 
to be published in the near future. This fiber was over 0.5 km long and 
was characterized by an 18-um-diameter core, a 15-um cladding thick- 
ness, a 100-um overall diameter, and a numerical aperture of 0.17. 
Loss minima occurred at 0.86, 0.90, and 1.02 um. The average losses 
at these wavelengths were 1.9, 2.4, and 1.1 dB/km, respectively. The 
loss at 1.06 um was 1.2 dB/km. 

In addition to low loss, optical waveguides should exhibit low pulse 
dispersion so that high data rates can be achieved. One way to accom- 
plish this is through the use of graded-refractive-index cores.’ By 
gradually changing the concentration of reactive gases as the film 
thickness is built up, graded-refractive-index profiles can be achieved. 
This has been accomplished in the GeO2-doped-core system by varying 
the concentration of GeCl, in the gas stream and in the silica-core, 
borosilicate-clad system by varying the concentration of BCl; during 
the deposition. The loss spectrum of a fiber in which the BO; concen- 
tration gradually changes from 0 at the center to about 20 percent at 
the core-cladding interface is presented in Fig. 3. This fiber had a 
22-um core diameter, a 15-um cladding thickness, and a 0.17 numerical 
aperture. Minima occur in the loss spectrum at 0.90 and 1.04 yum. 
The losses at these wavelengths were 2.3 and 1.7 dB/km, respectively. 
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